转录组入门(6)： reads计数
点击此处查看最新的网赚项目教程
要求
实现这个功能的软件也很多，还是烦请大家先自己搜索几个教程，入门请统一用htseq-count，对每个样本都会输出一个表达量文件。
需要用脚本合并所有的样本为表达矩阵。参考：生信编程直播第四题：多个同样的行列式文件合并起来
对这个表达矩阵可以自己简单在excel或者R里面摸索，求平均值，方差。
看看一些生物学意义特殊的基因表现如何，比如GAPDH,β-ACTIN等等。
理论基础
在上篇的比对中，我们需要纠结是否真的需要比对，如果你只需要知道已知基因的表达情况，那么可以选择alignment-free工具（例如salmon, sailfish)，如果你需要找到noval isoforms，那么就需要alignment-based工具（如HISAT2, STAR）。到了这一篇的基因（转录本）定量，需要考虑的因素就更加多了，以至于我不知道如何说清才能理清逻辑。
定量分为三个水平
在基因水平上，常用的软件为HTSeq-count，featureCounts，BEDTools, Qualimap, Rsubread, GenomicRanges等。以常用的HTSeq-count为例，这些工具要解决的问题就是根据read和基因位置的overlap判断这个read到底是谁家的孩子。值得注意的是不同工具对multimapping reads处理方式也是不同的，例如HTSeq-count就直接当它们不存在。而Qualimpa则是一人一份，平均分配。
对每个基因计数之后得到的count matrix再后续的分析中，要注意标准化的问题。如果你要比较同一个样本(within-sample)不同基因之间的表达情况，你就需要考虑到转录本长度，因为转录本越长，那么检测的片段也会更多，直接比较等于让小孩和大人进行赛跑。如果你是比较不同样本（across sample）同一个基因的表达情况，虽然不必在意转录本长度，但是你要考虑到测序深度（sequence depth)，毕竟测序深度越高，检测到的概率越大。除了这两个因素外，你还需要考虑GC%所导致的偏差，以及测序仪器的系统偏差。目前对read count标准化的算法有RPKM（SE）, FPKM（PE），TPM, TMM等，不同算法之间的差异与换算方法已经有文章进行整理和吐槽了。但是，有一些下游分析的软件会要求是输入的count matrix是原始数据，未经标准化，比如说DESeq2，这个时候你需要注意你上一步所用软件会不会进行标准化。
在转录本水平上，一般常用工具为Cufflinks和它的继任者StringTie， eXpress。这些软件要处理的难题就时转录本亚型（isoforms）之间通常是有重叠的，当二代测序读长低于转录本长度时，如何进行区分？这些工具大多采用的都是expectation maximization（EM）。好在我们有三代测序。
上述软件都是alignment-based，目前许多alignment-free软件，如kallisto, silfish, salmon，能够省去比对这一步，直接得到read count，在运行效率上更高。不过最近一篇文献[1]指出这类方法在估计丰度时存在样本特异性和读长偏差。
在外显子使用水平上，其实和基因水平的统计类似。但是值得注意的是为了更好的计数，我们需要提供无重叠的外显子区域的gtf文件[2]。用于分析差异外显子使用的DEXSeq提供了一个Python脚本（dexseq_prepare_annotation.py）执行这个任务。
小结
计数分为三个水平： gene-level, transcript-level, exon-usage-level
标准化方法： FPKM RPKM TMM TPM
输出表达矩阵
在RNA-Seq分析中，每一个基因就是一个feature（特征？），而基因被认为是它的所有外显子的和集。在可变剪切分析中，可以单独把每个外显子当作一个feature。而在ChIP-Seq分析中，feature则是预先定义的结合域。但是确定一个read到底属于哪一个feature有时会非常棘手。因此HTSeq提供了三种模式，示意图见前一幅图
基本用法非常的简单：
# 安装
conda install htseq
# 使用
# htseq-count [options]
htseq-count -r pos -f bam RNA-Seq/aligned/SRR3589957_sorted.bam reference/gencode.v26lift37.annotation.sorted.gtf > SRR3589957.count用一个循环处理多个BAM文件(在/mnt/f/Data目录下)
mkdir -p RNA-Seq/matrix/
for i in `seq 56 58`
do
htseq-count -s no -r pos -f bam RNA-Seq/aligned/SRR35899${i}_sorted.bam reference/gencode.v26lift37.annotation.sorted.gtf > RNA-Seq/matrix/SRR35899${i}.count 2> RNA-Seq/matrix/SRR35899${i}.log
done运行的时间会比较久，所以可以去了解不同参数的用法了，其中比较常用的为：
gene_id “ENSG00000223972.5_2”; transcript_id “ENST00000456328.2_1”; gene_type “transcribed_unprocessed_pseudogene”; gene_name “DDX11L1”; transcript_type “processed_transcript”; transcript_name “DDX11L1-002”; exon_number 2; exon_id “ENSE00003582793.1_1”; level 2;…
Jimmy的伏笔
我们这次分析是人类mRNA-Seq测序的结果，但是我们其实只下载了3个sra文件。一般而言RNA-Seq数据分析都至少要有2个重复，所以必须要有4个sra文件才行。我在仔细读完文章的方法这一段以后，发现他们有一批数据用的是其他课题组的： For 293 cells, the mRNA-seq results of the control samples include (1) those from the doxycycline-treated parental Flp-In T-REx 293 cells by us and (2) those from the doxycycline-treated control Flp-In T-REx 293 cells performed by another group unrelated to us (sample GSM1095127 in GSE44976)。 然后和Jimmy交流之后，他也承认自己只分析了小鼠的数据，而没有分析人类的数据。所以我们需要根据文章提供的线索下载另外一份数据，才能进行下一步的分析。
这个时候就有一个经常被问到的问题：不同来源的RNA-Seq数据能够直接比较吗？甚至说如果不同来源的RNA-seq数据的构建文库都不一样该如何比较?不同来源的RNA-Seq结果之间比较需要考虑 批次效应（batch effect) 的影响。
处理批次效应，根据我搜索的结果，是不能使用FPKM/RPKM，关于这个标准化的吐槽，我在biostars上找到了如下观点：
有人建议使用一个Bioconductor包 我没有具体了解，有生之年去了解补充。
还有人引用了一篇文献 IVT-seq reveals extreme bias in RNA-sequencing 证明不同文库的RNA-Seq结果会存在很大差异。
结论： 可以问下原作者他们是如何处理数据的，居然有一个居然没有重复的分析也能过审。改用小鼠数据进行分析。或者使用无重复的分析方法，或者模拟一份数据出来，先把流程走完。
合并表达矩阵
HTSeq-count输出结果是一个个独立的文件，后续分析需要把多个文件合并成一个行为基因名，列为样本名，中间为count的行列式文件。肯定是不会用Excel手动处理（虽然可以写一个VBA脚本，但是数据量过大不好处理了），这里使用的Python写一个脚本。
基本逻辑：
读取文件
建立一个字典，如果key不在字典中，新增key和value，如果key在字典中，追加value。
输出
#!/usr/bin/python3
import sys
mydict = {}
for file in sys.argv[1:]:
for line in open(file, 'r'):
key,value = line.strip().split('t')
if key in mydict:
mydict[key] = mydict[key] + 't' + value
else:
mydict[key] = value
for key,value in mydict.items():
print(key + 't' + value +'n')这几行代码写了2个番茄钟，但是debug花了我一个番茄钟。问题出在str和list两种数据格式的混乱使用。还有一个bug： 由于词典是无序的，所以原本代表样本来源的第一行，会跑到其他行。
在论坛上找到一个更加简洁的代码（要求基因名顺序排列），用到paste, awk, printf这三个shell命令。
paste *.txt | awk '{printf $1 "t";for(i=2;i<=NF;i+=2) printf $i"t";printf $i}'保存为countCombiner.py，输入文件为count, 输出为标准输出，需要重定向。
简单分析
这一步需要用到R语言或者是Excel读取数据。
1.导入数据
options(stringsAsFactors = FALSE)# import data if sample are smallcontrol <- read.table("F:/Data/RNA-Seq/matrix/SRR3589956.count",
sep="t", col.names = c("gene_id","control"))
rep1 <- read.table("F:/Data/RNA-Seq/matrix/SRR3589957.count",
sep="t", col.names = c("gene_id","rep1"))
rep2 <- read.table("F:/Data/RNA-Seq/matrix/SRR3589958.count",
sep="t",col.names = c("gene_id","rep2"))2.数据整合。gencode的注释文件中的gene_id（如ENSG00000105298.13_3）在EBI是不能搜索到的，所以我就只保留ENSG00000105298这部分。
# merge data and delete the unuseful row
raw_count <- merge(merge(control, rep1, by="gene_id"), rep2, by="gene_id")
raw_count_filt <- raw_count[-1:-5,]
ENSEMBL <- gsub("(.*?)\.\d*?_\d", "\1", raw_count_filt$gene_id)
row.names(raw_count_filt) <- ENSEMBL3.总体情况, 大部分基因都为0，所以可以删掉节省体积。
summary(raw_count_filt)
control rep1 rep2
Min. : 0.0 Min. : 0.0 Min. : 0.0
1st Qu.: 0.0 1st Qu.: 0.0 1st Qu.: 0.0
Median : 0.0 Median : 0.0 Median : 0.0
Mean : 356.1 Mean : 370.3 Mean : 316.6
3rd Qu.: 15.0 3rd Qu.: 15.0 3rd Qu.: 10.0
Max. :161867.0 Max. :121902.0 Max. :105565.04.看几个具体基因
在EBI上搜索GAPDH找到ID为ENSG00000111640。
GAPDH <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000111640",]
gene_id control rep1 rep2
ENSG00000111640 ENSG00000111640.14_2 41857 53902 55302文章研究的AKAP95（ENSG00000105127）的表达量在KD中都降低了
> AKAP95 AKAP95
gene_id control rep1 rep2
ENSG00000105127 ENSG00000105127.8_2 1168 539 506下面的差异基因表达，让我想想，该如何收拾Jimmy挖的坑。
参考文献
[1] Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis
[2] Detecting differential usage of exons from RNA-seq data
[3] RNA-seq Data Analysis-A Practical Approach(2015)
———END———
限 时 特 惠： 本站每日持续更新海量各大内部创业教程，一年会员只需98元，全站资源免费下载 点击查看详情
站 长 微 信： qs62318888
主题授权提示：请在后台主题设置-主题授权-激活主题的正版授权，授权购买：RiTheme官网
