芯片探针序列的基因组注释
这是我第二次在标题上写重磅！价值一千元的代码，虽然下面的技能或者说代码对我来说是非常简单啦，但是在有需求的粉丝看来真正的价值不可估量。
第一次是：
纯粹的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技巧啦。
■■ ■
最新文章
-
王长绅紫微斗数之廉贞天相在子午守财帛宫，命宫为武曲星
19分钟前 -
磁共振扫描操作磁共振冠脉成像（一）
33分钟前 -
“病急乱投医”
50分钟前 -
择日——吊用“太阳、太阴”吉神到宫的方法
1小时前 -
古诗词中的十大人生境界，古今一付笑谈中
1小时前 -
史上最全K线图大全：70种+42个K线组合图解
2小时前 -
明朝最大遗憾：崇祯若答应李自成条件，明朝不会灭亡，还灭了满清
2小时前 -
孔子语录100句（附释文）国学经典
2小时前
热门文章
-
业医必修——孙思邈《大医精诚》
1周前 -
怎样从门口看房屋吉凶
1周前 -
低估了广东的“甘蔗马蹄茅根水”，喝完半个小时就舒服了！
1周前 -
痰湿体质如何调理 最全痰湿体质调理方法(8)
1周前 -
自制家庭小咸菜
1周前 -
七彩神龙指标源码（八戒神龙主图无未来函数）
7天前 -
《奇门遁甲预测术》奇门：格局庚
6天前 -
10种朝鲜咸菜的腌制方法
7天前
有话要说...