# 下載gtf：選擇的是ensembl的GRCh38.94.gtf：
下載地址：ftp://ftp.ensembl.org/pub/release-94/gtf/homo_sapiens/
下載命令：
axel -n 10 ftp://ftp.ensembl.org/pub/release-94/gtf/homo_sapiens/Homo_sapiens.GRCh38.94.gtf.gz
# 解壓：
gunzip Homo_sapiens.GRCh38.94.gtf.gz
# head查看一下gtf的格式：
# 提取exon出來：
awk '{if ($3=="exon") print $0}' Homo_sapiens.GRCh38.94.gtf? > exon.gtf
# 去掉雙引號：
sed -i 's/"http://g' exon.gtf
# 提取需要的信息到txt：
awk 'BEGIN{FS="\t| |;";OFS="\t"}{print $1,$3,$4,$5,$7,$10,$25,$5-$4+1}' exon.gtf > exon.txt
其實這里提取genename/geneid、genelength這兩項就夠了。
# 查看下有多少exon，多少個gene：
wc -l exon.txt
-> 1262162
awk '{print $6}' exon.txt | sort | uniq |wc -l
-> 58735 (gene_id)
awk '{print $7}' exon.txt | sort | uniq |wc -l
-> 57169 (gene_name)
看來gene_id的數(shù)量比genename要多，用name吧。
# 計算長度：
awk 'BEGIN{OFS="\t"}{ s[$7] += $8 }END{for (i in s){print i,s[i]}} ' exon.txt > gene.length.txt
wc -l gene.length.txt
-> 57169 (gene_name)
同樣的思路也可以對轉(zhuǎn)錄本id/name的長度進行統(tǒng)計。
