以前做基因注释总用R包:TxDb.Hsapiens.UCSC.hg38.knownGene,但是今天做发现1个问题,就是有的基因特别长,与ncbi中检索出来的不一样,查找原因,在生物技能树上也看到了相关的问题,于是是果断参考用了gencode.v27.annotation.gtf.gz。
之前做法:
#有问题
Gene=as.data.frame(genes(TxDb.Hsapiens.UCSC.hg38.knownGene))
ann <- tryCatch(
suppressWarnings(select(org.Hs.eg.db,
keys=unique(Gene$gene_id),
keytype="ENTREZID",
columns=c("SYMBOL"))))
colnames(ann)[1]="gene_id"
Gene=merge(Gene,ann)
gn=as(Gene,"GRanges")
************************************************
现在参考基因组文件及注释:
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_27/gencode.v27.annotation.gtf.gz
grep -v "^#" gencode.v27.annotation.gtf|awk '$3=="gene"{print $1,$2,$3,$4,$5,$7,$12,$14}'|sed 's/;/\t/g'|sed 's/"//g'|sed 's/ /\t/g' >gencode.v2.annotation.gene
Gene=read.table("gencode.v2.annotation.gene",sep="\t",header=F,stringsAsFactors = F) gn=GRanges(seqnames =Gene$V1,ranges = IRanges(start=Gene[,4],end=Gene$V5),strand = Gene$V6,Symbols=Gene$V9)
peak1=GRanges(seqnames = paste("chr",sv_dt$Chr1,sep=""),ranges=IRanges(sv_dt$Pos1,sv_dt$Pos1),mcols=sv_dt[,1:2]) over.Peak1=mergeByOverlaps(peak1,gn)