2019-08-28 基因注释

以前做基因注释总用R包：TxDb.Hsapiens.UCSC.hg38.knownGene，但是今天做发现１个问题，就是有的基因特别长，与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)