shell玩转各个不同的gtf和gff文件--ID转换杂谈

玩转gtf

查看gtf每一行的信息

zcat Homo_sapiens.GRCh38.86.gtf.gz |head -100|grep -v ^# |cut -f 2
zcat Homo_sapiens.GRCh38.86.gtf.gz |head -100|grep -v ^# |cut -f 9

发现第九行是基因ID

zcat Homo_sapiens.GRCh38.86.gtf.gz |head -100|grep -v ^# |cut -f 9|cut -d ';' -f 1 > gene_id.txt

第九行后用';'分隔后的第一行是gene_id

image.png

下载refseq gff文件

https://www.ncbi.nlm.nih.gov/projects/genome/guide/human/index.shtml

zcat GRCh38_latest_genomic.gff.gz |head -100|grep -v '^#'|cut -f 1,9

第一列和第九列有gene_id信息

image.png

zcat GRCh38_latest_genomic.gff.gz |head -100|grep -v '^#'|cut -f 9|cut -d ';' -f 5
zcat GRCh38_latest_genomic.gff.gz |head -100|grep -v '^#'|cut -f 9|cut -d ';' -f 5|grep 'LOC'

第九列的按‘；’分隔后的第五列包含了基因名信息

image.png

还是没达到我要的目的，所以我下载了genecode的注释信息

参考资料https://zhuanlan.zhihu.com/p/36275161
gtf下载链接https://www.gencodegenes.org/human/
ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_30/gencode.v30.annotation.gtf.gz

zcat gencode.v30.annotation.gtf.gz |head -100|grep -v '^#'|cut -f 9

同样是第九列包含了

提取gencode内的基因信息

step1 探索数据

zcat gencode.v30.annotation.gtf.gz |head -100|grep -v '^#'| cut -f 9| cut -d ';' -f 1 |cut -d '"' -f 2

image.png

step2 paste - - 将两个信息贴一起

zcat gencode.v30.annotation.gtf.gz |grep -v '^#'| cut -f 9| cut -d ';' -f 4|cut -d '"' -f2 > gene_name
zcat gencode.v30.annotation.gtf.gz |head -100|grep -v '^#'| cut -f 9| cut -d ';' -f 1 |cut -d '"' -f 2 > ENSG
 paste ENSG gene_name > ID_tran

image.png

step3 将level 去掉并去重

cat ID_tran |head -100| grep -v 'level'
cat ID_tran | grep -v 'level'| uniq -c > ID_trans.txt

最后回到R里操作

rm(list = ls())
options(stringsAsFactors = F)
expr<- read.table(file = './our_result.csv.txt', sep = '\t', header = T)
ID_trans<- read.table(file = '../../gtf and gff/ID_trans.txt')
save(expr,ID_trans,file = 'step0.Rdata')


load(file = 'step0.Rdata')
table(expr$Gene%in%ID_trans$V3)
library(stringr)
ID_trans$gene<- str_split(ID_trans$V3,pattern = '-',simplify = T)[,1]
table(expr$Gene%in%ID_trans$gene)
ID_trans$gene<- str_split(ID_trans$V3,pattern = '[.]',simplify = T)[,1]
table(expr$Gene%in%ID_trans$gene)
table(expr$Gene%in%ID_trans$V2)

image.png