转录组入门(6): reads计数

要求

实现这个功能的软件也很多,还是烦请大家先自己搜索几个教程,入门请统一用htseq-count,对每个样本都会输出一个表达量文件。
需要用脚本合并所有的样本为表达矩阵。参考:生信编程直播第四题:多个同样的行列式文件合并起来
对这个表达矩阵可以自己简单在excel或者R里面摸索,求平均值,方差。
看看一些生物学意义特殊的基因表现如何,比如GAPDH,β-ACTIN等等。

理论基础

在上篇的比对中,我们需要纠结是否真的需要比对,如果你只需要知道已知基因的表达情况,那么可以选择alignment-free工具(例如salmon, sailfish),如果你需要找到noval isoforms,那么就需要alignment-based工具(如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则是一人一份,平均分配。

_images/count_modes.png

对每个基因计数之后得到的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提供了三种模式,示意图见前一幅图

  • the union of all the sets S(i) for mode union. This mode is recommended for most use cases.
  • the intersection of all the sets S(i) for mode intersection-strict.
  • the intersection of all non-empty sets S(i) for mode intersection-nonempty.

基本用法非常的简单:

# 安装
conda install htseq
# 使用
# htseq-count [options] <alignment_file> <gtf_file>
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

运行的时间会比较久,所以可以去了解不同参数的用法了,其中比较常用的为:

  • -f bam/sam: 指定输入文件格式,默认SAM
  • -r name/pos: 你需要利用samtool sort对数据根据read name或者位置进行排序,默认是name
  • -s yes/no/reverse: 数据是否来自于strand-specific assay。DNA是双链的,所以需要判断到底来自于哪条链。如果选择了no, 那么每一条read都会跟正义链和反义链进行比较。默认的yes对于双端测序表示第一个read都在同一个链上,第二个read则在另一条链上。
  • -a 最低质量, 剔除低于阈值的read
  • -m 模式 union(默认), intersection-strict and intersection-nonempty。一般而言就用默认的,作者也是这样认为的。
  • -i id attribute: 在GTF文件的最后一栏里,会有这个基因的多个命名方式(如下), RNA-Seq数据分析常用的是gene_id, 当然你可以写一个脚本替换成其他命名方式。

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上找到了如下观点:

  • 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 used because 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 because gene length is constant

有人建议使用一个Bioconductor包http://www.bioconductor.org/packages/devel/bioc/html/sva.html 我没有具体了解,有生之年去了解补充。

还有人引用了一篇文献 IVT-seq reveals extreme bias in RNA-sequencing 证明不同文库的RNA-Seq结果会存在很大差异。

结论: 可以问下原作者他们是如何处理数据的,居然有一个居然没有重复的分析也能过审。改用小鼠数据进行分析。或者使用无重复的分析方法,或者模拟一份数据出来,先把流程走完。

合并表达矩阵

HTSeq-count输出结果是一个个独立的文件,后续分析需要把多个文件合并成一个行为基因名,列为样本名,中间为count的行列式文件。肯定是不会用Excel手动处理(虽然可以写一个VBA脚本,但是数据量过大不好处理了),这里使用的Python写一个脚本。

基本逻辑:

  1. 读取文件
  2. 建立一个字典,如果key不在字典中,新增key和value,如果key在字典中,追加value。
  3. 输出
#!/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 small
control <- 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) <- 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-seq analysis

[2] Detecting differential usage of exons from RNA-seq data

[3] RNA-seq Data Analysis-A Practical Approach(2015)

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 148,637评论 1 318
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 63,443评论 1 266
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 99,164评论 0 218
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 42,075评论 0 188
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 50,080评论 1 266
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 39,365评论 1 184
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 30,901评论 2 283
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 29,649评论 0 176
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 33,122评论 0 223
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 29,734评论 2 225
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,093评论 1 236
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 27,548评论 2 222
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,028评论 3 216
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,765评论 0 9
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,291评论 0 178
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 34,162评论 2 239
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 34,293评论 2 242

推荐阅读更多精彩内容