搜索
查看: 3427|回复: 1

[mRNA-seq] 转录组作业(六)reads计数

[复制链接]

11

主题

22

帖子

330

积分

中级会员

Rank: 3Rank: 3

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

基因(转录本定量)
定量分为三个水平:基因水平(gene-level)转录本水平(transcript-level)外显子使用水平(exon-usage-level)
基因水平
基因水平上,常用的软件为HTSeq-count,featureCounts,BEDTools, Qualimap, Rsubread, GenomicRanges等,主要用来解决的问题是根据read和基因位置的overlap判断这个read到底是谁家的孩子,不同的工具对multimapping reads处理方式也是不同的,例如 HTSeq-count直接当他们不存在。而  Qualimap则是一人一份,平均分配。
对每个基因计数之后得到的count matrix 再后续的分析中,要注意标准化的问题。
比较同一样本(within-sample)不同基因之间的表达情况,需要考虑到转录本长度
比较不同样本(across sample)同一基因表达情况,不必在意转录本长度,但是要考虑到测序深度
除此之外,还必须要考虑到GC%所导致的偏差,以及测序仪器的系统偏差。read count标准化的算法有RPKM(SE), FPKM(PE),TPM, TMM。有一些下游分析的软件会要求是输入的count matrix 是原始数据,未经过标准化,比如说DESeq2,这个时候需要注意上一步使用软件会不会进行标准化。
步骤:累计比对到原始reads数(HTSeq),基因水平的定量(使用GTF注释文件),基因表达量归一算法(RPKM,FPKM,TPM)
RPKM:Reads Per Kilobase per Million mapped reads的缩写,RPKM是将map到基因的read数除以map到基因组上的所有read数(以million为单位)与RNA的长度(以KB为单位)。RNA-seq是二代测序技术中用来表示基因表达量或丰度的方法。在衡量基因表达量时,若是单纯以map到的read数来计算基因的表达量,在统计上是不合理的。因为在随机抽样的情况下,序列较长的基因被抽到的机率本来就会比序列短的基因较高,如此一来,序列长的基因永远会被认为表达量较高,而错估基因真正的表现量  。FPKM:Fragments Per Kilobase of transcript per Million fragments mapped,每1百万个map上的reads中map到map到外显子的每1k个碱基上的fragments个数。
FPKM和RPKM计算方法基本抑制,主要一个采用reads一个采用Fragments,片段比读段含义更广。因此FPKM包含含义更广。
RPKM不仅对测序深度做了归一化,也对基因长度做了归一化,使得不同长度的基因在不同测序深度下得到的基因表达水平估计具有了可比性,是目前最为常用的基因估计方法。
TPM:   Transcripts per million 每百万读段中来自某基因的读段数,考虑了测序深度对读段计数的影响,但没有基因组长度的影响,主要用于miRNA测序。
转录本水平
转录本水平上,  一般常用工具为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)执行这个任务。
Q:外显子使用水平主要用来做什么呢?
输出表达矩阵
在RNA-Seq分析中,每一个基因就是一个feature(特征?),而基因被认为是它的所有外显子的和集。在可变剪切分析中,可以单独把每个外显子当作一个feature。(而在ChIP-Seq分析中,feature则是预先定义的结合域。)CHIP-Seq没有接触过。但是确定一个read到底属于哪一个feature有时会非常棘手。因此HTSeq提供了三种模式,示意图见前一幅图。
处理多个bam文件 for i in `seq 56 58`
do
htseq-count -s no -r pos -f bam SRR35899${i}_sorted.bam /mnt/d/RNA_seq/data/reference/genome/hg19true/gencode.v26lift37.annotation.sorted.gtf > matrix/SRR35899${i}.count 2> matrix/SRR35899${i}.log
done
-s yes/no/reverse:数据是否来自strand-specific assay,选择no:每一条read都会跟正义链和反义链进行比较,默认yes对于双端测序表示第一个read都在同一条链上,第二个read则在另一条链上。
-r name/pos:默认name,利用samtool sort 对数据根据name还是position进行排序。(Panda姐说用name比较好,因为曾经他用pos出现错误)
-f bam/sam:指定输入文件格式,默认SAM
-a 最低质量, 剔除低于阈值的read
-m 模式 union(默认), intersection-strict and intersection-nonempty。一般而言就用默认的,作者也是这样认为的。
-i id attribute: 在GTF文件的最后一栏里,会有这个基因的多个命名方式(如下), RNA-Seq数据分析常用的是gene_id, 当然你可以写一个脚本替换成其他命名方式。
小插曲:出现错误:  Error occurred when processing SAM input (record #1 in file /mnt/d/AKAP95/aligned/SRR3589959.nsorted.bam):   my_showwarning() takes from 3 to 5 positional arguments but 6 were given   [Exception type: TypeError, raised in warnings.py:99]  
解决方案:是因为htseq与python版本不和导致的,详见http://www.biotrainee.com:8080/thread-2086-1-1.html       #激活python=2.7的环境
conda info --envs
source activate sunshine
#下载安装
conda install -c bioconda htseq=0.7.2
批次效应
涉及问题:不同来源的RNA-Seq数据能够直接比较吗?甚至说如果不同来源的RNA-seq数据的构建文库都不一样该如何比较?不同来源的RNA-Seq结果之间比较需要考虑 批次效应(batch effect) 的影响。这一步我也很关心,在芯片处理中可以用R包来处理批次效应。
合并表达矩阵
   用R和perl
perl -lne 'if($ARGV=~/SRR35899(.*).count/){print "$1\t$_"}' *| grep -v EnsEMBL_Gene_ID > hg.count
#在R中,把工作目录设为文件所在
hg <- read.csv(file = "hg.count",header = F,sep = "\t")
colnames(hg) <- c('sample','gene','count')
library(reshape2)
reads <- dcast(hg,formula = gene ~ sample)
write.table(reads,file = "hg_join.count",sep = "\t",quote = FALSE,row.names = FALSE)出错的地方1
代码如下samtools view -S SRR3589956.sam -b > SRR3589956.bam
samtools sort SRR3589956.bam -o SRR3589956_sorted.bam
samtools index SRR3589956_sorted.bam
#reads计数
htseq-count -s no -f bam SRR3589956_sorted.bam /mnt/d/RNA_seq/data/reference/genome/hg19true/gencode.v26lift37.annotation.sorted.gtf > matrix/SRR3589956.count
#输出Warning: 60287703 reads with missing mate encountered.
62939247 SAM alignment pairs processed.
htseq版本:0.7.2检查:Total records:                          64852965                 总比对记录有64852965QC failed:                              0                  
Optical/PCR duplicate:                  0                         PCR重复数
Non primary hits                        7139405                 多位点匹配
Unmapped reads:                         2628387             不匹配的read数
mapq < mapq_cut (non-unique):           4173091     mapq >= mapq_cut (unique):              50912082
Read-1:                                 25654550
Read-2:                                 25257532
Reads map to '+':                       25493062               比对到+链
Reads map to '-':                       25419020                比对到-链
Non-splice reads:                       40663134               没有剪切位点的reads数
Splice reads:                           10248948                  有剪接位点的reads数
Reads mapped in proper pairs:           49490102      匹配的reads数
Proper-paired reads map to different chrom:0总体匹配率还是有70%以上的啊。查了网上的,说是可以忽略。。。不可能吧。后面又加载了最新版本的htseq,一样的结果。htseq加了-r pos 以后出现此问题,求助orz Mate pairing was ambiguous for 46484 records; mate key for first such record: ('SRR3589958.sra.6645700', 'first', 'chr1', 16914, 'chr1', 17008, 145).
29843679 SAM alignment pairs processed.
出错的地方2
reads的名字就是不标准。所以导致后面导入的R里面的时候,检索不出来。嗯,待解决。




本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:转录组作业(5)序列比对
下一篇:R circos
回复

使用道具 举报

11

主题

22

帖子

330

积分

中级会员

Rank: 3Rank: 3

积分
330
 楼主| 发表于 2017-11-6 17:08:50 | 显示全部楼层
后面我仔细看了大家的作业,也都是存在这个问题的。
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-8-22 03:04 , Processed in 0.031046 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.