问题描述链特异文库鉴定长链非编码RNA(lncRNA)的基本步骤是
- stringtie每个样本分别组装得到转录本，获得是一个gtf文件
stringtie -p 12 --rf -o transcripts.gtf sorted.bam
ls *.gtf > list.merged.txt
~/Biotools/gffcompare/gffcompare -r reference.gtf -i list.merged.txt -o merged
得到一个 merged.combined.gtf这个文件里给每一个转录本分配了一个class_code用来表示转录本相对于参考基因组的位置
以上图片来源于论文 GFF Utilities: GffRead and GffCompare
长链非编码RNA通常是选择class_code为u/x/i，比如论文 Global identification of Arabidopsis lncRNAs revealsthe regulation of MAF4 by a natural antisense RNA中就提到 only transcripts with TAIR10 annotation [Cufflinks class codes ‘u’ (intergenic transcripts),’x’ (Exonic overlap with reference on the opposite strand),’i’ (transcripts entirely within intron) were retained.
那么问题就来了，如何利用 merged.combined.gtf这个文件获得 class_code 为 u、x和i的转录本的gtf文件呢
找到了一个办法，python中有一个模块 pyGTF，github链接是https://github.com/chengcz/pyGTF
直接使用pip安装
pip install pyGTF
可以解析gft格式的注释文件
利用这个模块来写一个简单的脚本
import sys
from pyGTF import Transcript
from pyGTF import GTFReader
in_gft = sys.argv[1]
class_code = sys.argv[2]
out_gtf = sys.argv[3]
fw = open(out_gtf,'w')
with GTFReader(in_gtf,flag_stream=True) as fi:
for i in fi:
if i._attri['class_code'] == class_code:
i.to_gtf(fw)
fw.close()
使用方法是
python 01.py in.gtf i out.gtf
####今天学到的另外一个知识点： samtools统计fasta文件序列长度，根据序列名提取序列
参考https://www.cnblogs.com/xudongliang/p/5200655.html
使用命令
samtools faidx input.fasta
会生成一个input.fasta.fai的文件，文件的内容总共有5列 第一列是序列名，第二列是序列长度，第四列是每行多少个碱基
根据序列名提取序列 这里好像只能提取单条序列
samtools faidx input.fasta TCONS_00000018 > TCONS_00000018.fa
还可以加上指定的位置
samtools faidx input.fasta TCONS_00000018:1-10
>TCONS_00000018:1-10
TGGGCGAACG
