开发者社区> 徐洲更> 正文
阿里云
为了无法计算的价值
打开APP
阿里云APP内打开

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

简介: 要求 实现这个功能的软件也很多,还是烦请大家先自己搜索几个教程,入门请统一用htseq-count,对每个样本都会输出一个表达量文件。 需要用脚本合并所有的样本为表达矩阵。
+关注继续查看

要求

实现这个功能的软件也很多,还是烦请大家先自己搜索几个教程,入门请统一用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则是一人一份,平均分配。

img_57c7e5330aed0ad76b6cab783db20260.jpe
_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)

版权声明:本文内容由阿里云实名注册用户自发贡献,版权归原作者所有,阿里云开发者社区不拥有其著作权,亦不承担相应法律责任。具体规则请查看《阿里云开发者社区用户服务协议》和《阿里云开发者社区知识产权保护指引》。如果您发现本社区中有涉嫌抄袭的内容,填写侵权投诉表单进行举报,一经查实,本社区将立刻删除涉嫌侵权内容。

相关文章
LIO-SAM代码逐行解读(4)-IMU预积分(1)
LIO-SAM代码逐行解读(4)-IMU预积分(1)
0 0
LIO-SAM代码逐行解读(4)-IMU预积分(2)
LIO-SAM代码逐行解读(4)-IMU预积分(2)
0 0
基于分数FLMS 和可变功率分数 LMS (VP-FLMS) 实现系统识别附matlab代码
基于分数FLMS 和可变功率分数 LMS (VP-FLMS) 实现系统识别附matlab代码
0 0
Seurat 4.0 | 单细胞转录组数据整合(scRNA-seq integration)
Seurat 4.0 | 单细胞转录组数据整合(scRNA-seq integration)
0 0
如何使用GMAP/GSNAP进行转录组序列比对
GMAP最早用于讲EST/cDNA序列比对到参考基因组上,可以用于基因组结构注释。后来高通量测序时代,又开发了GSNAP支持高通量数据比对,这篇文章主要介绍GMAP,毕竟高通量转录组数据比对大家更喜欢用STAR, HISTA2等软件。
1631 0
用fastp对转录组数据做QC
链接 fastp: 极速全能的FASTQ文件自动质控+过滤+校正+预处理软件 github地址 看到介绍的时候是真的心动不已↓↓↓ fastp可以仅仅扫描 FASTQ 文件一次,就完成比FASTQC + cutadapt + Trimmomatic 这...
1311 0
转录组入门(5): 序列比对
欢迎来GitHub上fork,一起进步: https://github.com/xuzhougeng 比对软件很多,首先大家去收集一下,因为我们是带大家入门,请统一用hisat2,并且搞懂它的用法。
1700 0
转录组入门(2):读文章拿到测序数据
本系列课程学习的文章是:AKAP95 regulates splicing through scaffolding RNAs and RNA processing factors.
834 0
基因表达分析(前传)-准备count矩阵
还在利用hisat, tophat这些耳熟能详的软件将read比对到基因组(转录组)上,然后统计每个基因的count数么?试试这些不需要比对,速度更快的工具吧。
1420 0
+关注
徐洲更
生信媛公众号编辑、生信必修课之软件安装课程作者
文章
问答
文章排行榜
最热
最新
相关电子书
更多
低代码开发师(初级)实战教程
立即下载
阿里巴巴DevOps 最佳实践手册
立即下载
冬季实战营第三期:MySQL数据库进阶实战
立即下载