芯片探針序列的基因組注釋
這是我第二次在标題上寫重磅！價值一千元的代碼，雖然下面的技能或者說代碼對我來說是非常簡單啦，但是在有需求的粉絲看來真正的價值不可估量。
第一次是：
純粹的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技巧啦。
■■ ■
最新文章
-
李辛最新采訪 | 孩子身心問題背後的原因（下）
12分鐘前 -
針家心悟：針灸如何取穴（純幹貨）
32分鐘前 -
有哪些非常實用的職場道理，是當了領導以後才明白的
50分鐘前 -
世界地理分區之歐洲篇
1小時前 -
陳益峰：飛鵝地，王氏正二品武将的老祖墳
1小時前 -
狗子從小頂白菜，頂草莓，結果頂着頂着頂出一個冠軍
2小時前 -
直覺主義——數學概念是自主的智力活動
2小時前 -
張本美和欲重金聘請張繼科執教
2小時前
熱門文章
-
威士忌高階關于 OB 與 IB，一次給你說透！
1天前 -
為什麼五點鐘要起床答案讓人吃驚！（現在知道還不晚）
1周前 -
美麗中國-2870：中國最大的内陸河，塔裡木河
3天前 -
老張老李侃門球之140篇
3天前 -
這個穴位可以治療多種胃痛腹痛，還可以減肥
1周前 -
診餘雜記（師傳經驗）
5天前 -
二十四山開門放水作竈真訣開門放水作竈直訣——子山
3天前 -
門球技巧 隻需五個字讓你打好門球
3天前
有話要說...