用gffread从基因组中提取基因、转录本、蛋白、启动子、非编码

https://jishuin.proginn.com/p/763bfbd5e36e
1.获取转录本序列
gffread GRCh38.gtf -g GRCh38.fa -w GRCh38.transcripts.fa

2. 获取CDS序列
   gffread GRCh38.gtf -g GRCh38.fa -x GRCh38.cds.fa
3. 获取蛋白序列
   gffread GRCh38.gtf -g GRCh38.fa -y GRCh38.protein.fa
   4.提取基因启动子序列
   首先确定启动子区域，这里定义转录起始位点上游1000 bp和下游500 bp为启动子区域。
   sed 's/"/\t/g' GRCh38.gtf | awk 'BEGIN{OFS=FS="\t"}{if(7=="+") {start=4+500;} else {if(5-500; end=1,start,end,10,$7;}}' >GRCh38.promoter.bed

然后提取序列。这里用到了bedtools工具，官方有提供编译好的二进制文件，下载下来即可使用。
-name: 输出基因名字（bed文件的第四列）
-s: 考虑到正反链（对于启动子区域，是否考虑链的信息关系不太大）
bedtools getfasta -name -s -fi GRCh38.fa -bed GRCh38.promoter.bed >GRCh38.promoter.fa
5.提取基因序列
提取基因序列的操作也类似于提取启动子序列。这里要注意GFF文件的序列位置是从1开始，而bed文件的位置是从0开始，前闭后开，所以要对序列的起始位置进行-1的操作。
type="gene"
sed 's/"/\t/g' GRCh38.gtf | awk -v type="3==type) {print 4-1,14,".",{type}" 'BEGIN{OFS=FS="\t"}{if(1,5,20,3=="exon") {tr=5; oldtr=tr;} else {print 4-1,tr,7; start=$5;} } }' >GRCh38.intron.bed