怎么根據(jù)class_code篩選轉(zhuǎn)錄本

鏈特異文庫鑒定長鏈非編碼RNA(lncRNA)的基本步驟是

• hisat2將原始測序數(shù)據(jù)比對到參考基因組
• samtools獲得排序的bam文件
• stringtie每個樣本分別組裝得到轉(zhuǎn)錄本，獲得是一個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這個文件里給每一個轉(zhuǎn)錄本分配了一個class_code用來表示轉(zhuǎn)錄本相對于參考基因組的位置

以上圖片來源于論文 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的轉(zhuǎn)錄本的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 

####今天學(xué)到的另外一個知識點： samtools統(tǒng)計fasta文件序列長度，根據(jù)序列名提取序列

https://www.cnblogs.com/xudongliang/p/5200655.html

使用命令

samtools faidx input.fasta

會生成一個input.fasta.fai的文件，文件的內(nèi)容總共有5列 第一列是序列名，第二列是序列長度，第四列是每行多少個堿基

根據(jù)序列名提取序列 這里好像只能提取單條序列

samtools faidx input.fasta TCONS_00000018 > TCONS_00000018.fa

還可以加上指定的位置

samtools faidx input.fasta TCONS_00000018:1-10
>TCONS_00000018:1-10
TGGGCGAACG