芯片探针序列的基因组注释

这是我第二次在标题上写重磅！价值一千元的代码，虽然下面的技能或者说代码对我来说是非常简单啦，但是在有需求的粉丝看来真正的价值不可估量。

第一次是：

纯粹的R代码技巧，怕粉丝看不懂，我已经花了一个星期做铺垫：

• 1

• 2

• 3

• 4

• 5

• 6

前面我提到过有些芯片，各种地方都是找不到其设计的探针对应的基因的，但是探针序列一般会给出，比如：

HumanLncRNAExpressionArrayV4.0AS-LNC-H-V4.020,730mRNAsand40,173LncRNAs8*60K

以前我会简单的回答，其实就是芯片探针的重新注释，重点是

• probe sequences 探针序列下载

• uniquely mapped to the human genome (hg19) by Bowtie without mismatch. 参考基因组下载及比对

• chromosomal position of lncRNA genes based on annotations from GENCODE (Release 23)坐标提取，最后使用bedtools进行坐标映射

三部曲罢了，不过感觉会linux的朋友不多，我还是用R来一波这个操作。

首先下载序列

这里我选择在GEO官网的GPL平台下载 : https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL21827

rm(list=ls())##魔幻操作，一键清空~
options(stringsAsFactors=F)
#注意查看下载文件的大小，检查数据
f='GPL21827_eSet.Rdata'

library(GEOquery)
#这个包需要注意两个配置，一般来说自动化的配置是足够的。
#Settingoptions('download.file.method.GEOquery'='auto')
#Settingoptions('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
gset<-getGEO('GPL21827',destdir=".")##平台文件
save(gset,file=f)##保存到本地
}
load('GPL21827_eSet.Rdata')##载入数据
class(gset)
length(gset)
gset
colnames(Table(gset))
probe2symbol=Table(gset)[,c(1,4)]
all_recs=paste(apply(probe2symbol,1,function(x)paste0('>',x[1],'\n',x[2])),collapse='\n')
temp<-tempfile()##编程技巧，把变量写入临时文件~
write(all_recs,temp)

这个技巧我在生信菜鸟团博客发布过：http://www.bio-info-trainee.com/3732.html 芯片概况如下：

然后对人类的参考基因组构建索引并且比对

主要是参考基因组下载会耗费时间，还有构建索引耗时也很可观！

library(Rsubread)
#推荐从ENSEMBL上面下载成套的参考基因组fa及基因注释GTF文件
dir='~/data/project/qiang/release1/Genomes/'
ref<-file.path(dir,'Homo_sapiens.GRCh38.dna.toplevel.fa')
buildindex(base,reference=ref)
##是单端数据，fa序列来源于上一个步骤输出的gpl的探针
reads<-temp
align(index="reference_index",readfile1=reads,
output_file="alignResults.BAM",phredOffset=64)
propmapped("alignResults.BAM")

构建大约耗时一个小时，具体如下：

比对速度很快，因为探针序列只有6万多，如下：

在linux下得到比对后的bam文件也很简单的。

读入人类基因组注释文件

也是需要一点点R技巧，可以参考我在生信菜鸟团的博客：http://www.bio-info-trainee.com/3742.html 使用refGenome加上dplyr玩转gtf文件

library(Rsubread)
#推荐从ENSEMBL上面下载成套的参考基因组fa及基因注释GTF文件
dir='~/data/project/release1/Genomes/'
gtf<-file.path(dir,'Homo_sapiens.GRCh38.82.gtf')
if(!require(refGenome))install.packages("refGenome")
#createensemblGenomeobjectforstoringEnsemblgenomicannotationdata
ens<-ensemblGenome()
#readGTFfileintoensemblGenomeobject
read.gtf(ens,useBasedir=F,gtf)

class(ens)
#countsallannotationsoneachseqname
tableSeqids(ens)
#returnsallannotationsonmitochondria
extractSeqids(ens,'MT')
#summarisefeaturesinGTFfile
tableFeatures(ens)
#createtableofgenes
my_gene<-getGenePositions(ens)
dim(my_gene)

#lengthofgenes
gt=my_gene
my_gene_length<-gt$end-gt$start
my_density<-density(my_gene_length)
plot(my_density,main='Distributionofgenelengths')
##重点是要成为对象
library(GenomicRanges)
my_gr<-with(my_gene,GRanges(seqid,IRanges(start,end),
strand,id=gene_id))

如果是linux的shell脚本，一句话就可以搞定其实。

坐标映射

把自己制作好的bam文件的坐标，跟提取自gtf文件的坐标信息对应起来，使用GenomicRanges包自带的函数即可。

值得注意的是把bam文件读入R，并且转为grange对象需要一点点技巧，我在生信菜鸟团博客写过：http://www.bio-info-trainee.com/3740.html

library(Rsamtools)
bamFile="alignResults.BAM"
quickBamFlagSummary(bamFile)
#https://kasperdanielhansen.github.io/genbioconductor/html/Rsamtools.html
bam<-scanBam(bamFile)
bam
names(bam[[1]])
tmp=as.data.frame(do.call(cbind,lapply(bam[[1]],as.character)))
tmp=tmp[tmp$flag!=4,]#60885probes
#intersect()ontwoGRangesobjects.
library(GenomicRanges)
my_seq<-with(tmp,GRanges(as.character(rname),
IRanges(as.numeric(pos)-60,as.numeric(pos)+60),
as.character(strand),
id=as.character(qname)))
gr3=intersect(my_seq,my_gr)
gr3
o=findOverlaps(my_seq,my_gr)
o
lo=cbind(as.data.frame(my_seq[queryHits(o)]),
as.data.frame(my_gr[subjectHits(o)]))
head(lo)
write.table(lo,file='GPL21827_probe2ensemb.csv',row.names=F,sep=',')

当然，坐标映射本身也是满满的R技巧啦。

■■ ■