转录组⼊门（6）：reads计数
要求
小学英语学习网站实现这个功能的软件也很多，还是烦请⼤家先⾃⼰搜索⼏个教程，⼊门请统⼀⽤htq-count，对每个样本都会输出⼀个表达量⽂件。
需要⽤脚本合并所有的样本为表达矩阵。参考：⽣信编程直播第四题：多个同样的⾏列式⽂件合并起来
reveals
对这个表达矩阵可以⾃⼰简单在excel或者R⾥⾯摸索，求平均值，⽅差。
看看⼀些⽣物学意义特殊的基因表现如何，⽐如GAPDH,β-ACTIN等等。
理论基础
在上篇的⽐对中，我们需要纠结是否真的需要⽐对，如果你只需要知道已知基因的表达情况，那么可以选择alignment-free⼯具（例如salmon, sailfish)，如果你需要找到noval isoforms，那么就需要alignment-bad⼯具（如HISAT2, STAR）。到了这⼀篇的基因（转录本）定量，需要考虑的因素就更加多了，以⾄于我不知道如何说清才能理清逻辑。
定量分为三个⽔平
基因⽔平(gene-level)
转录本⽔平(transcript-level)
外显⼦使⽤⽔平(exon-usage-level)。
在基因⽔平上，常⽤的软件为HTSeq-count，featureCounts，BEDTools, Qualimap, Rsubread, GenomicRanges等。以常⽤的HTSeq-count为例，这些⼯具要解决的问题就是根据read和基因位置的overlap判断这个read到底是谁家的孩⼦。值得注意的是不同⼯具对multimapping reads处理⽅式也是不同的，例如HTSeq-count就直接当它们不存在。⽽Qualimpa则是⼀⼈⼀份，平均分配。cyclobenzaprine
_images/count_modes.png
对每个基因计数之后得到的count matrix再后续的分析中，要注意标准化的问题。如果你要⽐较同⼀个样本(within-sample)不同基因之间
的表达情况，你就需要考虑到转录本长度，因为转录本越长，那么检测的⽚段也会更多，直接⽐较等于让⼩孩和⼤⼈进⾏赛跑。如果你是⽐
较不同样本（across sample）同⼀个基因的表达情况，虽然不必在意转录本长度，但是你要考虑到测序深度（quence depth)，毕竟
名侦探柯南口头禅
测序深度越⾼，检测到的概率越⼤。除了这两个因素外，你还需要考虑GC%所导致的偏差，以及测序仪器的系统偏差。⽬前对read count
标准化的算法有RPKM（SE）, FPKM（PE），TPM, TMM等，不同算法之间的差异与换算⽅法已经有⽂章进⾏整理和吐槽了。但是，有
⼀些下游分析的软件会要求是输⼊的count matrix是原始数据，未经标准化，⽐如说DESeq2，这个时候你需要注意你上⼀步所⽤软件会不
会进⾏标准化。
在转录本⽔平上，⼀般常⽤⼯具为Cufflinks和它的继任者StringTie， eXpress。这些软件要处理的难题就时转录本亚型（isoforms）之
间通常是有重叠的，当⼆代测序读长低于转录本长度时，如何进⾏区分？这些⼯具⼤多采⽤的都是expectation maximization（EM）。好
在我们有三代测序。
上述软件都是alignment-bad，⽬前许多alignment-free软件，如kallisto, silfish, salmon，能够省去⽐对这⼀步，直接得到read
count，在运⾏效率上更⾼。不过最近⼀篇⽂献[1]指出这类⽅法在估计丰度时存在样本特异性和读长偏差。
在外显⼦使⽤⽔平上，其实和基因⽔平的统计类似。但是值得注意的是为了更好的计数，我们需要提供⽆重叠的外显⼦区域的gtf⽂件[2]。
⽤于分析差异外显⼦使⽤的DEXSeq提供了⼀个Python脚本（dexq_prepare_annotation.py）执⾏这个任务。
⼩结
计数分为三个⽔平： gene-level, transcript-level, exon-usage-level
标准化⽅法： FPKM RPKM TMM TPM
输出表达矩阵
在RNA-Seq分析中，每⼀个基因就是⼀个feature（特征？），⽽基因被认为是它的所有外显⼦的和集。在可变剪切分析中，可以单独把每dhc是什么
个外显⼦当作⼀个feature。⽽在ChIP-Seq分析中，feature则是预先定义的结合域。但是确定⼀个read到底属于哪⼀个feature有时会⾮
常棘⼿。因此HTSeq提供了三种模式，⽰意图见前⼀幅图
the union of all the ts S(i) for mode union. This mode is recommended for most u cas.
the interction of all the ts S(i) for mode interction-strict.
the interction of all non-empty ts S(i) for mode interction-nonempty.
基本⽤法⾮常的简单：
# 安装
conda install htq
# 使⽤
# htq-count [options] <alignment_file> <gtf_file>
htq-count -r pos -f bam RNA-Seq/aligned/SRR3589957_sorted.bam reference/gencode.v26lift37.f > unt
⽤⼀个循环处理多个BAM⽂件(在/mnt/f/Data⽬录下)
mkdir -p RNA-Seq/matrix/
for i in `q 56 58`
do
htq-count -s no -r pos -f bam RNA-Seq/aligned/SRR35899${i}_sorted.bam reference/gencode.v26lift37.f > RNA-Seq/matrix/SRR35899${ done
运⾏的时间会⽐较久，所以可以去了解不同参数的⽤法了，其中⽐较常⽤的为：
最适合你的英文名字-f bam/sam： 指定输⼊⽂件格式，默认SAM
-r name/pos: 你需要利⽤samtool sort对数据根据read name或者位置进⾏排序，默认是name
-s yes/no/rever: 数据是否来⾃于strand-specific assay。DNA是双链的，所以需要判断到底来⾃于哪条链。如果选择了no， 那么每⼀条read都会跟正义链和反义链进⾏⽐较。默认的yes对于双端测序表⽰第⼀个read都在同⼀个链上，第⼆个read则在另⼀条链
上。
-a 最低质量， 剔除低于阈值的read
-m 模式 union（默认）, interction-strict and interction-nonempty。⼀般⽽⾔就⽤默认的，作者也是这样认为的。
-i id attribute: 在GTF⽂件的最后⼀栏⾥，会有这个基因的多个命名⽅式（如下）， RNA-Seq数据分析常⽤的是gene_id， 当然你可以写⼀个脚本替换成其他命名⽅式。
gene_id "ENSG00000223972.5_2"; transcript_id "ENST00000456328.2_1"; gene_type
"transcribed_unprocesd_pudogene"; gene_name "DDX11L1"; transcript_type "procesd_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-q results of the control samples include (1) tho from the doxycycline-treated parental Flp-In T-REx 293 cells by us and (2) tho 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-q数据的构建⽂库都不⼀样该如何⽐较?不同来源的RNA-Seq结果之间⽐较需要考虑 批次效应（batch effect) 的影响。
处理批次效应，根据我搜索的结果，是不能使⽤FPKM/RPKM，关于这个标准化的吐槽，我在biostars上找到了如下观点：FPKM/RPKM 不是标准化的⽅法，它会引⼊⽂库特异的协变量
FPKM/RPKM has never been peer-reviewed, it has been introduced as an ad-hoc measure in a supplementary 没有同⾏评审
One of the authors of this paper states, that it should not be ud becau of faulty arithmetic 作者说算法有问题
All reviews so far have shown it to be an inferior scale for DE analysis of genes Length normalization is mostly
dispensable imo in DE analysis becau gene length is constant
还有⼈引⽤了⼀篇⽂献 IVT-q reveals extreme bias in RNA-quencing 证明不同⽂库的RNA-Seq结果会存在很⼤差异。
结论： 可以问下原作者他们是如何处理数据的，居然有⼀个居然没有重复的分析也能过审。改⽤⼩⿏数据进⾏分析。或者使⽤⽆重复的分析⽅法，或者模拟⼀份数据出来，先把流程⾛完。
合并表达矩阵
HTSeq-count输出结果是⼀个个独⽴的⽂件，后续分析需要把多个⽂件合并成⼀个⾏为基因名，列为样本名，中间为count的⾏列式⽂件。肯定是不会⽤Excel⼿动处理（虽然可以写⼀个VBA脚本，但是数据量过⼤不好处理了），这⾥使⽤的Python写⼀个脚本。
基本逻辑：
1. 读取⽂件
2. 建⽴⼀个字典，如果key不在字典中，新增key和value，如果key在字典中，追加value。
3. 输出
#!/usr/bin/python3
import sys
mydict = {}
boulder
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
el:
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 small
control <- read.table("F:/Data/RNA-Seq/unt",
p="\t", col.names = c("gene_id","control"))
rep1 <- read.table("F:/Data/RNA-Seq/unt",
vossp="\t", col.names = c("gene_id","rep1"))
rep2 <- read.table("F:/Data/RNA-Seq/unt",
p="\t",col.names = c("gene_id","rep2"))
2.数据整合。gencode的注释⽂件中的gene_id（如ENSG00000105298.13_3）在EBI是不能搜索到的，所以我就只保留
ENSG00000105298这部分。
# merge data and delete the unuful 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) <- ENSEMBL
3.总体情况, ⼤部分基因都为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.0
4.看⼏个具体基因
在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 <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000105127",]
> 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-q analysis
[2] Detecting differential usage of exons from RNA-q data
[3] RNA-q Data Analysis-A Practical Approach(2015)
