# 下载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/"//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的数量比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)
同样的思路也可以对转录本id/name的长度进行统计。