从NCBI上下载转录组数据,很多文章方法描述简单,无法判断是否为链特异性数据,导致在mapping和raw reads count时参数不知如何选择。一直弄不懂链特异性和参数设置的同志们可以去看《链特异建库那点事》,讲的非常清楚(虽然小白看完还是一知半解)。所以关于判断转录组数据是否为链特异性,偷懒的小白找到了一个比较省事儿的方法,用RSeQC的infer_experiment.py工具。
官网链接:http://rseqc.sourceforge.net/#infer-experiment-py
官网用法如下:需要输入自己数据的bam文件,这个容易拿到。
但是另一个hg19.refseq.bed12文件,官网给出的参数解释如下:但是Bedops(Bedopts)的gff2bed工具和RSeQC的举例bed12格式却大不一样。
Bedops(Bedopts)的gff2bed转化成bed文件的结果:
原链接:
https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/gtf2bed.html
而RSeQC举例的bed12格式的文件,除了包含常用的chromosome, start, end, name, score, strand等信息外,最后一列包含了多个extron和intron等的位置,用逗号隔开。
原链接:
http://dldcc-web.brc.bcm.edu/lilab/liguow/RSeQC/dat/sample.bed
小白通过搜索终于找到了获取bed12文件Reference gene model in bed format的方法,就是使用UCSC的gtfToGenePre工具,小白在上一篇笔记《Reference gene model in bed format》中已经详细讲述,这里只列代码:
#安装gtfToGenePre
conda install -c bioconda ucsc-gtftogenepred
#准备好基因组gtf文件,从gtf转换为GenePred格式
gtfToGenePred -genePredExt -geneNameAsName2 genes.gtf gene.tmp
#从GenePred文件提取信息得到bed文件
awk '{print $2"\t"$4"\t"$5"\t"$1"\t0\t"$3"\t"$6"\t"$7"\t0\t"$8"\t"$9"\t"$10}' gene.tmp > genes_refseq.bed12
拿到bed12文件后,开始试试用infer_experiment.py判断数据是否为链特异性。
infer_experiment.py -r genes_refseq.bed12 -i col.bam
结果:
Reading reference gene model genes_refseq.bed12 ... Done
Loading SAM/BAM file ... Total 200000 usable reads were sampled
This is SingleEnd Data
Fraction of reads failed to determine: 0.0050
Fraction of reads explained by "++,--": 0.4974
Fraction of reads explained by "+-,-+": 0.4976
这表明该数据为单端数据,以illumina standard建库方式为代表的fr-unstranded的非链特异性转录组数据。
再来一个链特异性的双端数据试试:
infer_experiment.py -r genes_refseq.bed12 -i m1_col_1_s.bam
结果:
Reading reference gene model genes_refseq.bed12 ... Done
Loading SAM/BAM file ... Total 200000 usable reads were sampled
This is PairEnd Data
Fraction of reads failed to determine: 0.0057
Fraction of reads explained by "1++,1--,2+-,2-+": 0.0049
Fraction of reads explained by "1+-,1-+,2++,2--": 0.9894
这表明该数据为双端数据,以dUTP建库方式为代表的RF (fr-firststrand)的链特异性转录组数据。
结果的详细解读可以去官网或参考链接:
http://rseqc.sourceforge.net/#infer-experiment-py
https://www.jianshu.com/p/4987dce4d165
欢迎关注和讨论哦,小白们一起学习生信~