有参转录组实战7-基因序列提取

#本教程仿自于“https://zhuanlan.zhihu.com/p/439168788”。

#正则表达式教程https://www.runoob.com/regexp/regexp-tutorial.html。

#1，提取转录本

gffread Ptri_genome.gtf -g Ptri_genome.fa -w Ptri.transcripts.fa

#2，CDS

gffread Ptri_genome.gtf -g Ptri_genome.fa -x Ptri.cds.fa

#3，Protein

gffread Ptri_genome.gtf -g Ptri_genome.fa -y Ptri.protein.fa

#4，Length of chromosomes

cut -f 1,2 Ptri_genome.fa.fai > Ptri_chr.size  

#5，检查GTF文件第一行

head -n 1 Ptri_genome.gtf | sed 's/"/\t/g' | tr '\t' '\n' | sed = | sed 'N;s/\n/\t/'

#6，检查GTF文件第二行

sed -n '2p' Ptri_genome.gtf | sed 's/"/\t/g' | tr '\t' '\n' | sed = | sed 'N;s/\n/\t/'

#7，提取启动子Promoter, 2000bp, first, create bed file

sed 's/"/\t/g' Ptri_genome.gtf | awk 'BEGIN{OFS=FS="\t"}{if($3=="transcript") {if($7=="+") {start=$4-2000; end=$4;} else {if($7=="-") start=$5; end=$5+2000; } if(start<0) start=0; print $1,start,end,$12,$12,$7;}}' >Ptri.promoter.bed

#8, Install bedtools

conda install bedtools

#9, Get the promoter sequences

bedtools getfasta -name -s -fi Ptri_genome.fa -bed Ptri.promoter.bed > Ptri.promoter.fa

#10, Simplify title but not recommend

cut -d ':' -f 1 Ptri.promoter.fa> Ptri.promoter.simplename.fa

#11, Gene, the bed file can view the gene location

type="transcript"

sed 's/"/\t/g' Ptri_genome.gtf | awk -v type="${type}" 'BEGIN{OFS=FS="\t"}{if($3==type) {print $1,$4-1,$5,$12,".",$7}}' > Ptri_gene.bed

#12, Get the gene sequences

bedtools getfasta -name -s -fi Ptri_genome.fa -bed Ptri_gene.bed > Ptri_gene.gene.fa

#13, Simplify title but not recommend

cut -d ':' -f 1 Ptri_gene.gene.fa > Ptri.gene.simplename.fa

#14, 最后保留这些文件，一些细节可用Notepad++ software修改，excel自行整合。不想理解代码，建议用TBtools，这软件在提取序列上还是很好用的。

#Ptri.transcripts.fa#转录本

#Ptri.cds.fa#CDS

#Ptri.protein.fa#蛋白

#Ptri_chr.size#染色体长度

#Ptri.promoter.fa#启动子

#Ptri_gene.bed#基因的位置信息

#Ptri_gene.gene.fa#基因序列

#虫师