PANDA姐的转录组入门(6)：reads计数

任务

• 前期补充：bam文件排序及sam/bam格式

• 搜索实现reads计数的软件教程

• 用htseq-count，每个样本都会输出一个表达量文件

• 编写脚本合并所有的样本为表达矩阵

• 对表达矩阵可进行简单的探索

• 看看一些生物学意义特殊的基因表现如何，比如GAPDH,β-ACTIN等等

代码记录

前期补充：bam文件排序及sam/bam格式

（要透彻地理解，一定要自己动手，光看不练是假把式！）

samtools sort 默认参数 
samtools view Homo_sapiens_AKAP95_KD_miR_12_293_cell_sorted.bam |less  -S 

看上去很类似fastq文件，它也有read名称，序列，质量等信息，但是又不完全一样。

首先，每个read只占一行，只是它被tab分成了很多列，一共有12列，分别记录了1：

1. read名称

2. SAM标记

3. chromosome

4. 5′端起始位置

5. MAPQ（mapping quality，描述比对的质量，数字越大，特异性越高）

6. CIGAR字串，记录插入，删除，错配以及splice junctions（后剪切拼接的接头）

7. mate名称，记录mate pair信息

8. mate的位置

9. 模板的长度

10. read序列

11. read质量

12. 程序用标记

显然，其中chromosome至CIGAR的信息都是非常重要的。但是这些对我们不重要，我们只需要了解SAM/BAM文件是什么，就可以了。重要的是如果进行下游的操作。更多详情见了解sam格式比对结果

1. samtools sort -n --threads 30 Homo_sapiens_AKAP95_KD_miR_12_293_cell.bam -o Homo_sapiens_AKAP95_KD_miR_12_293_cell_sorted_name.bam

2. samtools view Homo_sapiens_AKAP95_KD_miR_12_293_cell_sorted_name.bam |less  -S

1. # 按read name重新进行排序，后续reads count会更快，出错率低

2. ls --color=never  *.bam|while read id;do(samtools sort -n --threads 40  $id  -o ${id%.*}_sorted_name.bam);done

3. ls --color=never  *_sorted_name.bam | while read id;do(samtools index $id);done

4. # 不能根据reads name 排序的进行index

5. # [E::hts_idx_push] unsorted positions

6. # samtools index: "Homo_sapiens_AKAP95_KD_miR_12_293_cell_sorted_name.bam" is corrupted or unsorted

pos排序的文件小些，为什么呢？因为按照染色体位置排序相似的内容会排在一起，按照reads name排序内容则更加杂乱，跟为排序之前的内容差不多。所以一般没特殊需求，可以采用默认的pos排序。提高文件的压缩比例。

为什么 BAM 文件 sort 之后体积会变小 
因为BAM 文件是压缩的二进制文件，对文件内容排序之后相似的内容排在一起，使得文件压缩比提高了，因此排序之后的 BAM 文件变小了，相对应的 SAM 文件就是纯文本文件，对 SAM 文件进行排序就不会改变文件大小。而且由于 RNA-seq 中由于基因表达量的关系，RNA-seq 的数据比对结果 BAM 文件使用 samtools 进行 sort 之后文件压缩比例变化会比 DNA-seq 更甚。另外，samtools 对 BAM 文件进行排序之后那些没有比对上的 reads 会被放在文件的末尾。

搜索实现reads计数的软件教程

（终有一天，会派上大用场，用心学！）

From: Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis 

Figure 1 
The RNACocktail analysis protocol. RNACocktail is a comprehensive protocol of RNA-seq data analysis. The figure summarizes the widely used approaches for the key steps over the broad spectrum of RNA-seq analysis and also succinctly captures the possible workflows one can use to analyse RNA-seq data

• Salmon2

• bedtools multicov3

• HTseq4

1. usage: htseq-count [options] alignment_file gff_file

2. <alignment_file>：比对到基因组后得到的SAM文件 （SAMtools 包含一些perl脚本可以将大多数的比对文件转换成SAM格式 ），注意基因组mapping时一定要用支持剪接的比对软件（splicing-aware aligner）进行比对软件如TopHat.  HTSeq-count 需要用到SAM格式中的 CIGAR 区域的信息。

3. <gff_file>:  包含单位信息的gff/GTF文件（gff文件格式），大多数情况下就是指注释文件;   由于GTF文件其实就是gff文件格式的变形，在这里同样可以传入GTF格式文件。

4. This script takes one or more alignment files in SAM/BAM format and a feature file in GFF format and calculates for each feature the number of reads mapping to it. See http://htseq.readthedocs.io/en/master/count.html for details.

5. positional arguments:

6.  samfilenames          Path to the SAM/BAM files containing the mapped reads.

7.                        If '-' is selected, read from standard input 想要通过标准输入来传入基因组mapping得到SAM文件，用 – 替换 <alignment_file> 即可

8.  featuresfilename      Path to the file containing the features

9. optional arguments:

10.  -h, --help            show this help message and exit

11.  -f {sam,bam}, --format {sam,bam}

12.                        type of <alignment_file> data, either 'sam' or 'bam'

13.                        (default: sam)

14.                       指定输入文件的格式，<format> 可以是 sam (for text SAM files) 或 bam (for binary BAM files)格式，默认是sam.

15.  -r {pos,name}, --order {pos,name}

16.                        'pos' or 'name'. Sorting order of <alignment_file>

17.                        (default: name). **Paired-end sequencing data must be sorted either by position or by read name, and the sorting order must be specified. **Ignored for single-end data. 如果你是双端测序，必须要对SAM进行排序（单端可不必排序，但这里我也推荐对单端测序结果排序已减少内存消耗并提高软件效率），对read name或 染色体位置进行排序皆可（这里我推荐按read name排序，因为通过位置排序我遇到过错误）。具体需要通过 -r 参数指定，所以排序请详见参数 -r，在sort时用-n则不用修改，默认的排序则修改成pos。

18.  --max-reads-in-buffer MAX_BUFFER_SIZE

19.                        When <alignment_file> is paired end sorted by position, allow only so many reads to stay in memory until the mates are found (raising this number will use more memory). Has no effect for single end or paired end sorted by name

20.  -s {yes,no,reverse}, --stranded {yes,no,reverse}

21.                        whether the data is from a strand-specific assay.Specify 'yes', 'no', or 'reverse' (default: yes).

22.                        'reverse' means 'yes' with reversed strand interpretation

23.                         是否链特异性，如果不是修改成no

24.  -a MINAQUAL, --minaqual MINAQUAL

25.                        skip all reads with alignment quality lower than the given minimum value (default: 10)

26.                      指定一个最低 read mapping质量值，低于<minaqual>值会被过滤掉

27.  -t FEATURETYPE, --type FEATURETYPE

28.                        feature type (3rd column in GFF file) to be used, all features of other type are ignored (default, suitable for Ensembl GTF files: exon)

29.                        指定最小计数单位类型（gff文件的第3列中的类型如： exon）, 指定后其他单位类型将被忽略（默认情况下，对于rna-seq分析采用Ensembl GTF：http://mblab.wustl.edu/GTF22.html文件类型时,默认值是：exon）。

30.  -i IDATTR, --idattr IDATTR

31.                        GFF attribute to be used as feature ID (default,suitable for Ensembl GTF files: gene_id)

32.                       最终的计数单位，一般为基因。 默认为：gene_id，也可以设置转录本，但由于模型问题，计数效果不佳

33.                        需要改成gene_name，可以直接识别，如果是gene_id不一定认识

34.  --additional-attr ADDITIONAL_ATTR [ADDITIONAL_ATTR ...]

35.                        Additional feature attributes (default: none, suitable for Ensembl GTF files: gene_name)

36.  -m {union,intersection-strict,intersection-nonempty}, --mode {union,intersection-strict,intersection-nonempty}

37.                        mode to handle reads overlapping more than one feature

38.                        (choices: union, intersection-strict, intersection-nonempty; default: union)

39.                      判断一个reads属于某个基因的模型，用来判断统计reads的时候对一些比较特殊的reads定义是否计入。

40.                        <mode> 包括：默认的union和intersection-strict、 intersection-nonempty  （默认：union）。

41.  --nonunique {none,all}

42.                        Whether to score reads that are not uniquely aligned or ambiguously assigned to features

43.  -o SAMOUTS [SAMOUTS ...], --samout SAMOUTS [SAMOUTS ...]

44.                        write out all SAM alignment records into an output SAM file called SAMOUT, annotating each line with its feature assignment (as an optional field with tag 'XF')

45.                      输出所有alignment的reads到一个叫 <samout>的sam文件中，通过一个可选的sam标签 ‘ XF ’来标注每一行对应的单位和计数，可以不设置。

46.  -q, --quiet           suppress progress report 屏蔽程序报告和警告

47. Written by Simon Anders ([e-mail]), European Molecular Biology Laboratory (EMBL). (c) 2010. Released under the terms of the GNU General Public License v3. Part of the 'HTSeq' framework, version 0.8.0.

更详细参考htseq-count的用法5

如何判断一个reads属于某个基因， htseq-count 提供了union, intersection_strict, intersection_nonempty 3 种模型，如图（大多数情况下作者推荐用union模型）： 

如果这3种模型不能满足您的需求，您还可以通过htseq-count 自定义模型6。

用htseq-count，每个样本都会输出一个表达量文件

（参考别人的，仍要用心写自己的！）

使用版本：version 0.7.2

1. ls *.Nsort.bam |while read id;do (nohup samtools view $id | ~/.local/bin/htseq-count -f sam -s no -i gene_name - ~/reference/gtf/gencode/gencode.v25lift37.annotation.gtf 1>${id%%.*}.geneCounts 2>${id%%.*}.HTseq.log&);done

一个RNA-seq实战-超级简单-2小时搞定！7过程比较复杂，中间有个bam/sam转换过程，下面是我的命令8：

1. # 人类

2. # mkdir matrix

3. Homo_GTF=~/rna_seq/data/reference/genome/hg19/gencode.v26lift37.annotation.gtf

4. ls --color=never Homo_sapiens*sorted_name.bam|while read id;do(htseq-count  -f bam -i gene_name  -s no  $id  $Homo_GTF > matrix/${id%_*}.count  2> matrix/${id%_*}.log);done

5. # 老鼠

6. # 下载gtf：http://www.gencodegenes.org/mouse_stats/archive.html

7. # wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_mouse/release_M10/gencode.vM10.annotation.gtf.gz

8. # gzip -d gencode.vM10.annotation.gtf.gz

9. Mus_GTF=~/rna_seq/data/reference/genome/mm10/gencode.vM10.annotation.gtf

10. ls --color=never Mus_musculus*sorted_name.bam|while read id;do(htseq-count  -f bam -i gene_name  -s no  $id  $Mus_GTF > matrix/${id%_*}.count  2> matrix/${id%_*}.log);done

出现警告，我在网上检索了一下，应该没什么问题，可以忽略9 

编写脚本合并所有的样本为表达矩阵

（看，之前学的就用到了吧！）

wc -l Homo_sapiens*.count 
head -n 4 Homo_sapiens*.count 

具体参考生信编程直播第四题：多个同样的行列式文件合并起来10

1. perl -lne 'if($ARGV=~/Homo_sapiens_(.*)_sorted.count/){print "$1\t$_"}' *|grep -v  Homo_sapiens>hg.count

2. # 先把所有文件进行合并

3. setwd("~/rna_seq/work/matrix")

4. hg <- read.csv(file = "hg.count",header = F,sep = "\t")

5. colnames(hg) <- c('sample','gene','count')

6. library(reshape2)

7. reads <- dcast(hg,formula = gene ~ sample)

8. write.table(reads,file = "hg_join.count",sep = "\t",quote = FALSE,row.names = FALSE)

对表达矩阵可进行简单的探索

summary(reads) 

sample_mean <- group_by(hg,sample)%>% summarize(mean=mean(count))

gene_mean <- group_by(hg,gene)%>% summarize(mean=mean(count)) 

tmp <- full_join(hg,gene_mean,by='gene') 

看看一些生物学意义特殊的基因表现如何，比如GAPDH,β-ACTIN等等

（生物学意义，你懂么？）

---

1. #软件工具#比较Samtools和HTSeq读取和处理SAM/BAM格式文件 

2. 使用Bedtools对RNA-seq进行基因计数：http://www.bio-info-trainee.com/745.html 

3. 转录组HTseq对基因表达量进行计数：http://www.bio-info-trainee.com/244.html 

4. htseq-count的用法：http://yangl.net/2016/09/21/htseq-count_manual/ 

5. htseq官网：http://htseq.readthedocs.io/en/release_0.9.1/ 

6. 一个RNA-seq实战-超级简单-2小时搞定！：http://www.bio-info-trainee.com/2218.html 

7. htseq-counts跟bedtools的区别：http://www.bio-info-trainee.com/2022.html 

8. Question: Htseq Count:https://www.biostars.org/p/84907/ 

9. 生信编程直播第四题：多个同样的行列式文件合并起来：http://www.biotrainee.com/thread-603-1-1.html 

历史记录

~ 牛郎织女，千里相会呀 ~

<<  滑动查看下一张图片  >>