怎么根据class_code筛选转录本

链特异文库鉴定长链非编码RNA(lncRNA)的基本步骤是

• hisat2将原始测序数据比对到参考基因组
• samtools获得排序的bam文件
• stringtie每个样本分别组装得到转录本，获得是一个gtf文件

stringtie -p 12 --rf -o transcripts.gtf sorted.bam

• gffcompare合并所有样本的gtf文件

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