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

那个时候写教程，以软件安装，软件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)

相关推荐