生信入门参考资料

用RSeQC对比对后的转录组数据进行质控

2018-01-23  本文已影响109人  因地制宜的生信达人

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

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

详细列表如下：

数据库文件

RSeQC接受4种文件格式:

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

https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_GENCODE_v19_basic.Cancer_genes.bed.gz
https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_UCSC_knownGene.bed.gz
https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_RefSeq.bed.gz
https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_rRNA.bed.gz
https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_GENCODE_GENE_V19_comprehensive.bed.gz
https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19.HouseKeepingGenes.bed
https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_AceView.bed.gz
https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_Ensembl.bed.gz

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

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

软件安装

# 如果python版本没有问题，那么直接用pip即可安装
pip install RSeQC --user
# 如果有conda，那么更方便
conda install -c bioconda rseqc
## 依赖于python2.7
## 所以conda可能需要先创建python2.7的环境,再安装
conda info --envs
conda create -n py2.7 python=2.7 rseqc
source activate py2.7

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

cat bam_path.txt |while read id ; do
      echo $id
      file=$(basename $id )
      sample=${file%%.*}
      bam_stat.py -i $id >${sample}_bam_stat.log
      read_quality.py -i  $id -o ${sample}_read_quality 
      read_NVC.py -i  $id -o ${sample}_read_NVC 
      read_GC.py -i $id -o ${sample}_read_GC 
      read_duplication.py -i $id -o ${sample}_read_duplication 
done

#nohup  bash run.sh 1>run.log 2>&1 &
#source activate py2.7  
mkdir -p db
cd db
if [ ! -f hg19_RefSeq.bed ]; then
wget https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_RefSeq.bed.gz
wget https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg38_RefSeq.bed.gz
wget https://sourceforge.net/projects/rseqc/files/BED/Mouse_Mus_musculus/mm10_RefSeq.bed.gz

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

cd ../
bed='db/mm10_RefSeq.bed'
nohup geneBody_coverage.py -r $bed -i bam_path.txt  -o coverage 1>coverage.log 2>&1 & 
cat bam_path.txt |while read id ; do
      echo $id 
      file=$(basename $id )
      sample=${file%%.*} 
      junction_annotation.py -i   $id   -r $bed  -o  ${sample}_junction  
      RPKM_saturation.py -r $bed -d '1++,1--,2+-,2-+' -i $id -o ${sample}_RPKM_saturation
      read_distribution.py -i  $id  -r $bed  >${sample}_distribution.log
done

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

#==================================================
#All numbers are READ count
#==================================================

Total records:                          45722327

QC failed:                              0
Optical/PCR duplicate:                  0
Non primary hits                        2687735
Unmapped reads:                         2338796
mapq < mapq_cut (non-unique):           2045264

mapq >= mapq_cut (unique):              38650532
Read-1:                                 19631272
Read-2:                                 19019260
Reads map to '+':                       19320271
Reads map to '-':                       19330261
Non-splice reads:                       20690614
Splice reads:                           17959918
Reads mapped in proper pairs:           36737552
Proper-paired reads map to different chrom:0

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

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

Total Reads                   40695796
Total Tags                    64718115
Total Assigned Tags           61411678
=====================================================================
Group               Total_bases         Tag_count           Tags/Kb
CDS_Exons           34406515            45257520            1315.38
5'UTR_Exons         6859302             2274659             331.62
3'UTR_Exons         25952114            9778098             376.77
Introns             943281009           3254031             3.45
TSS_up_1kb          19391072            65573               3.38
TSS_up_5kb          88202906            155561              1.76
TSS_up_10kb         160360035           222457              1.39
TES_down_1kb        19659116            216878              11.03
TES_down_5kb        84349049            524626              6.22
TES_down_10kb       149723035           624913              4.17
=====================================================================

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

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

junction_annotation.py:

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

一般来说，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最明显).

上一篇下一篇

猜你喜欢

热点阅读