21-有参转录组实战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#基因序列

#虫师