生信分析、meta分析、数据挖掘
TCGA、GEO、SEER、Oncomine

还是用RSeQC对比对后的转录组数据做一下质控

早在3年前我就关注过这个软件,因为某个公司的转录组结题报告里面提到过,但那个时候我还是不是很懂NGS组学异同,不太理解质量控制,只知道流程,别人分析什么我就分析什么,所以那个时候我博客里面关于RSeQC的教程很有趣(http://www.bio-info-trainee.com/113.html) 复制粘贴该url到浏览器可以查看我当初的教程,或者点击最下面的阅读原文即可。

那个时候写教程,以软件安装,软件input和output为主,因为觉得新手最容易纠结的就是这些了,但是现在回过头来看,软件安装已经成了小菜一碟,对各种bam/sam/vcf/gtf也耳熟能详,所谓的input/output也不是问题了。

所以,再看看我最近是如何记录该软件的吧:

RSeQC包是一个python软件,最新版是 v2.6.4 , 依赖于: gcc; python2.7; numpy; R

它提供了一系列有用的小工具能够评估高通量测序尤其是RNA-seq数据,比如一些基本模块,检查序列质量, 核酸组分偏性, PCR偏性, GC含量偏性,还有RNA-seq特异性模块: 评估测序饱和度, 映射读数分布, 覆盖均匀性, 链特异性, 转录水平RNA完整性等

详细列表如下:

  • bam2fq.py

  • bam2wig.py

  • bam_stat.py

  • clipping_profile.py

  • deletion_profile.py

  • divide_bam.py

  • FPKM_count.py

  • geneBody_coverage.py

  • geneBody_coverage2.py

  • infer_experiment.py

  • inner_distance.py

  • insertion_profile.py

  • junction_annotation.py

  • junction_saturation.py

  • mismatchprofile.pynormalizebigwig.pyoverlay_bigwig.py

  • read_distribution.py

  • read_duplication.py

  • readGC.pyreadhexamer.py

  • read_NVC.py

  • read_quality.py

  • RNAfragmentsize.py

  • RPKMcount.pyRPKMsaturation.py

  • spilt_bam.py

  • splitpairedbam.py

  • tin.py

数据库文件

RSeQC接受4种文件格式:

  • BED 格式: Tab 分割, 12列的表示基因模型的纯文本文件

  • SAM 或BAM 格式: 用来存储reads 比对结果信息.

  • 染色体大小文件: 只有两列的纯文本文

  • Fasta文件的参考基因组

数据库文件根据参考基因组版本自行选择下载,我这里要下载的是hg19系列,下载地址如下:

  1. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_GENCODE_v19_basic.Cancer_genes.bed.gz

  2. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_UCSC_knownGene.bed.gz

  3. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_RefSeq.bed.gz

  4. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_rRNA.bed.gz

  5. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_GENCODE_GENE_V19_comprehensive.bed.gz

  6. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19.HouseKeepingGenes.bed

  7. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_AceView.bed.gz

  8. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_Ensembl.bed.gz

  9. ls *.gz|while read id ; do gunzip $id;done

希望读者能够明白,看教程一定要看规律,我为什么列出如此多的url,其实就是想你领悟它们的共性: https://sourceforge.net/projects/rseqc/files/ 你在浏览器打开就明白了。

### 软件安装

  1. # 如果python版本没有问题,那么直接用pip即可安装

  2. pip install RSeQC --user

  3. # 如果有conda,那么更方便

  4. conda install -c bioconda rseqc

  5. ## 依赖于python2.7

  6. ## 所以conda可能需要先创建python2.7的环境,再安装

  7. conda info --envs

  8. conda create -n py2.7 python=2.7 rseqc

  9. source activate py2.7

虽然该软件的使用命令非常多,但很多功能并不是用来诊断转录组测序的,所以不在我们的考虑范围内。下面是我们经常会用得到的:

  1. # nohup  bash run.sh 1>run.log 2>&1 &

  2. #source activate py2.7  

  3. mkdir -p db

  4. cd db

  5. if [ ! -f hg19_RefSeq.bed ]; then

  6. wget https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_RefSeq.bed.gz

  7. wget https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg38_RefSeq.bed.gz

  8. wget https://sourceforge.net/projects/rseqc/files/BED/Mouse_Mus_musculus/mm10_RefSeq.bed.gz

  9. ls *.gz|while read id ; do gunzip $id;done

  10. fi

  11. cd ../

  12. bed='db/mm10_RefSeq.bed'

  13. nohup geneBody_coverage.py -r $bed -i bam_path.txt  -o coverage 1>coverage.log 2>&1 &

  14. cat bam_path.txt |while read id ; do

  15.      echo $id

  16.      file=$(basename $id )

  17.      sample=${file%%.*}

  18.      junction_annotation.py -i   $id   -r $bed  -o  ${sample}_junction

  19.      bam_stat.py -i $id >${sample}_bam_stat.log  

  20.      RPKM_saturation.py -r $bed -d '1++,1--,2+-,2-+' -i $id -o ${sample}_RPKM_saturation

  21.      ## below just like fastqc

  22.      nohup read_quality.py -i  $id -o ${sample}_read_quality &

  23.      nohup read_NVC.py -i  $id -o ${sample}_read_NVC &

  24.      nohup read_GC.py -i $id -o ${sample}_read_GC &

  25.      nohup read_duplication.py -i $id -o ${sample}_read_duplication &

  26.      read_distribution.py -i  $id  -r $bed  >${sample}_distribution.log

  27. done  

用 bam_stat.py来统计总比对记录, PCR重复数Non Primary Hits表示多匹配位点, 不匹配的reads比对到+链的reads比对到-链的reads有剪切位点的reads等.

  1. #==================================================

  2. #All numbers are READ count

  3. #==================================================

  4. Total records:                          45722327

  5. QC failed:                              0

  6. Optical/PCR duplicate:                  0

  7. Non primary hits                        2687735

  8. Unmapped reads:                         2338796

  9. mapq < mapq_cut (non-unique):           2045264

  10. mapq >= mapq_cut (unique):              38650532

  11. Read-1:                                 19631272

  12. Read-2:                                 19019260

  13. Reads map to '+':                       19320271

  14. Reads map to '-':                       19330261

  15. Non-splice reads:                       20690614

  16. Splice reads:                           17959918

  17. Reads mapped in proper pairs:           36737552

  18. Proper-paired reads map to different chrom:0

可以看到比对效果非常赞,这个转录组很成功!

另外一个比较赞的小程序就是: read_duplication.py 结果一般如下:

  1. Total Reads                   40695796

  2. Total Tags                    64718115

  3. Total Assigned Tags           61411678

  4. =====================================================================

  5. Group               Total_bases         Tag_count           Tags/Kb

  6. CDS_Exons           34406515            45257520            1315.38

  7. 5'UTR_Exons         6859302             2274659             331.62

  8. 3'UTR_Exons         25952114            9778098             376.77

  9. Introns             943281009           3254031             3.45

  10. TSS_up_1kb          19391072            65573               3.38

  11. TSS_up_5kb          88202906            155561              1.76

  12. TSS_up_10kb         160360035           222457              1.39

  13. TES_down_1kb        19659116            216878              11.03

  14. TES_down_5kb        84349049            524626              6.22

  15. TES_down_10kb       149723035           624913              4.17

  16. =====================================================================

可以用一个饼图来表示,在生信技能树论坛里面还有人专门提问过。

用 geneBody_coverage.py来计算RNA-seq reads在基因上的覆盖度,这里推荐对所有的样本的 bam文件一起运行该程序进行诊断,如图:

junction_annotation.py:

输入一个 BAM或 SAM文件和一个 bed12格式的参考基因文件,这个模块将根据参考基因模型计算剪切融合(splice junctions)事件.

  • splice read: 一个RNA read,能够被剪切一次或多次

  • splice junction:多个跨越同一个内含子的剪切事件能够合并为一个 splicing junction.

一般来说,novel的junction区域总是有的,因为我们用的是ucsc的refseq参考注释集,本身就是不够完整的。

RPKM_saturation.py

  任何样本统计( RPKM)的精度受样本大小( 测序深度)的影响,重抽样切片是使用部分数据来评估样本统计量的精度的方法。这个模块从总的 RNA reads中重抽样并计算每次的 RPKM值,通过这样我们就能检测当前测序深度是不是够的(如果测序深度不够RPKM的值将不稳定,如果测序深度足够则RPKM值将稳定)。*默认情况下,这个模块将计算20个 RPKM值(分别是对个转录本使用5%,10%,…,95%的总 reads),所以非常消耗内存哦。


  在结果图中,Y轴表示 Percent Relative Error 或 Percent Error

说明:Q1,Q2,Q3,Q4是按照转录本表达量4分位分开的.Q1表示的是表达量低于25%的转录本,以此类推.可以看出:随着样本量升高, RPKM与实际值的偏差也在降低.而且转录本表达量越高这种趋势越明显(Q4最明显).

写在最后:

NGS组学分析流程的每一个步骤都应该是有充分的质量控制,主要是考虑到各个项目的实际情况可能会比较特殊,如果都走一样的自动化,流水线的流程,肯定是会有问题的。

明天给大家看看,问题主要是什么,敬请期待哈。

赞(0) 打赏
未经允许不得转载:医学SCI科研之家 » 还是用RSeQC对比对后的转录组数据做一下质控
分享到: 更多 (0)

评论 抢沙发

  • 昵称 (必填)
  • 邮箱 (必填)
  • 网址

meta分析、生信分析

meta、生信交流群综合科研交流群