我說(shuō)真的，在我寫(xiě)這篇以前我從來(lái)沒(méi)有想過(guò)「提取轉(zhuǎn)錄本」這種東西有什么好寫(xiě)的
直到我遇到了百合，這個(gè)破事也能卡我一整天
1 失敗的嘗試
岷江百合（L. regale）的1號(hào)染色體能有超過(guò)4GB的大小，這就意味著超過(guò)了2^31-1（約2.1Gb）
如果用gffread提取，會(huì)生成錯(cuò)誤的fai文件（染色體坐標(biāo)出現(xiàn)負(fù)數(shù)），因?yàn)橐绯隽?2bit的整形
即使人工替換正確的fai文件，也迷之原因無(wú)法正確定位
bedtools也迷之原因失敗
2 救星：AGAT
AGAT（Another Gff Analysis Toolkit）
這是一個(gè)perl寫(xiě)的能夠處理大基因組的工具
安裝AGAT
conda install -c bioconda agat
提取轉(zhuǎn)錄本
agat_sp_extract_sequences.pl --gff "/media/desk16/tl5024/Lregale/lrv2.gff" --fasta "/media/desk16/tl5024/Lregale/lrv2_mc.fa" -t mRNA -o mRNA.fasta --keep_attributes "ID"
如果要提取其他類(lèi)型的，換成cds或者gene、exon即可
提取出來(lái)的序列抬頭有一點(diǎn)復(fù)雜
>LregChr01G000010.mRNA1 gene=LregChr01G000010 seq_id=Chr01 type=cds
ATGGCCAGTCGCCGGCGAGCCCTCCTCAAAATCATCGTCCTCGGCGACAGCGGAGTCGGA
AAGACATCATTGCTCAATCAATATGTGCACAAGAAGTTTAGGCAGCAGTACAAAGCTACA
ATCGGCACAGATTTCCTCACCAAGGAGCTCCAGATGGAGGACAGGCTCGTGACGTTGCAA
ATATGGGACACAGCAGGCCAGGAAAGATTTCAGAGTCTCGGTGTGGCGTTCTACCGTGGG
GCGGATTGCTGTATTCTGGCCTATGATGTCAACGTGCGCAAGTCATTTGATACCCTTGAC
AACTGGCATGAAGAGTTTCTTAAACAGGCTAGGGTTTCTGATCCCAAGACTTTTCCATTT
GTTCTACTTGGAAACAAAGTTGACATAGATGGTGGCAACAGCAGAGTGACCTCGGAGAAG
AAAGCTGGAGAATGGTGTGCTTCCAAGGGCAGCATACCTTACTTTGAAACTTCAGCAAAA
GAGGATTATAATGTTGACGCTGCTTTCTTATGTGCTGCTAAACTTGCACTAGCGAAAGAT
因此需要對(duì)抬頭進(jìn)行一點(diǎn)精簡(jiǎn)
sed 's/^>\([^ ]*\).*/>\1/' mRNA.fasta > mRNA_simple_header.fasta
# simplified file will be saved in mRNA_simple_header.fasta
搞定！
