目的
一个新的物种（没有参考基因组）用Trinity从头组装得到的转录本fasta文件，需要将其转换成gtf文件
思路
- 用http://www.bioinformatics.nl/courses/RNAseq/1c%20Assembly.pdf的思路，用StringTie得到gtf
- 用http://seqanswers.com/forums/showthread.php?t=44455的思路
bwa-mem of Trinity.fasta to reference genome
samtools -Sbh | bam2bed | bed12ToBed6 | bedToGenePred stdin stdout | genePredToGtf file stdin assembly.gtf
思路的实操记录
- 用
Trinity组装成fasta文件
- 用
- cd-hit去掉上述fasta文件的冗余序列
教程：①教程 | 如何用cd-hit去除冗余序列？ ②无参转录组序列组装 — Trinity
- cd-hit去掉上述fasta文件的冗余序列
-
- 利用
hisat2-build构建索引 | 教程： RNASEQ分析入门笔记3-HISAT2建立索引并比对序列
hisat2-build -p 10 test.fasta prefix - 利用
-
- 利用以上述组装得到的fasta文件作为参考序列，用
hisat2将原始的fastq文件比对到该fasta文件上，得到sam文件
hisat2 -x prefix -1 raw1.fq.gz -2 raw2.fq.gz -S test.sam &>log.txt - 利用以上述组装得到的fasta文件作为参考序列，用
-
- 将sam文件转换成bam文件
nohup samtools view -bS test.sam > test.bam & -
- 将
bam格式转换成bed6格式
nohup bamToBed -i ./test.bam > test.bed & # 其实我到了这步已经是bed6文件格式了 nohup bed12ToBed6 -i test.bed > test.bed & # 将bed12转换成bed6 - 将
-
- 将
bed6格式转换成genePred格式
bedToGenePred ./test.bed ./test.genePred - 将
-
- 将
genePred格式转换成gtf格式
genePredToGtf file test.genePred test.gtf - 将
或者直接用管道写成一条命令也行
