内容不是教程，均为个人学习笔记，以备后日查询
学习资源均来自生信技能树论坛，生信技能树B站视频，微信公众号：生信技能树
为啥是我最不想用的方法？
- 常用平台可通过bioconductor中R包完成，详见生信技能树教程：用R获取芯片探针与基因的对应关系
- 没有对应的R包，也可以用Jimmy老师写的R包ANNOPROBE来完成，其中包含了大多数平台的探针信息。
以上两种方法都可以轻易的获取平台探针对应的gene_symbol。
但最悲伤的是，费劲儿学了R语言和一些生信分析技能能后，发现自己想分析的芯片对应的平台通过以上方法都无法获得探针信息，比如我要做的：GPL5175。连里面的基因是啥都不知道，就别分析了=。=！
只好跟着生信技能树的教程自己下载探针序列进行比对，电脑不行，相当费时。
- 16G内存的电脑跑下面的流程真的让我等的花儿都谢了。
- 有时间赶紧给linux学了，瘸腿跑的感觉很Der。
1.下载分析过程中需要用的文件
-
探针对应的序列
一般用getGEO()函数下载对应的平台文件就能获取探针与序列的对应数据。可惜，通过GEO数据库下载GPL5175中没有探针序列，只好去芯片的官网下载了，地址一般在GEO数据库中可见，如下：
-
参考基因组fa文件
一定要下Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz里面有个primary，别看错了。 -
基因注释GTF文件
两个都下载好：
2.将下载的好的探针序列整理成需要的格式，怎么整理就是R语言操作的事儿。格式如下：
3.人参考基因组构建索引并且比对
library(Rsubread)
dir='D:/Chrom download'' #这里就是Homo_sapiens.GRCh38.dna.chromosome.1.fa文件的所在路径
ref <- file.path(dir,'Homo_sapiens.GRCh38.dna.chromosome.1.fa')
buildindex(basename="reference_index",reference=ref)
reads <- temp
align(index="reference_index",readfile1=reads,
output_file="alignResults.BAM",phredOffset=64)
propmapped("alignResults.BAM")
无尽的等待，16G内存的忧伤=。=！
4.将基因组注释文件载入R
library(Rsubread)
dir='D:/Chrom download' #这里就是Homo_sapiens.GRCh38.99.gtf文件的所在路径
gtf <- file.path(dir,'Homo_sapiens.GRCh38.99.gtf')
library(refGenome)
ens <- ensemblGenome()
read.gtf(ens,useBasedir = F,gtf)
extractSeqids(ens, 'MT')
tableFeatures(ens)
my_gene <- getGenePositions(ens)
gt<-my_gene
my_gene_length <- gt$end - gt$start
my_density <- density(my_gene_length)
library(GenomicRanges)
my_gr <- with(my_gene, GRanges(seqid, IRanges(start, end),
strand, id = gene_id))
这里有一点需要注意，refGenome这个R包已经从CRAN repository移除了，通过install()函数是无法安装的，需要手动下载压缩包然后安装。
5.坐标映射
library(Rsamtools)
bamFile="alignResults.BAM"
quickBamFlagSummary(bamFile)
bam <- scanBam(bamFile)
tmp=as.data.frame(do.call(cbind,lapply(bam[[1]], as.character)))
tmp=tmp[tmp$flag!=4,]
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)
o = findOverlaps(my_seq,my_gr)
lo=cbind(as.data.frame(my_seq[queryHits(o)]),
as.data.frame(my_gr[subjectHits(o)]))
write.table(lo,file = 'GPL5175_probe2ensemb.csv',row.names = F,sep = ',')#保存结果
再使用R读进来
GPL5175_probe2ensemb<-read.csv('GPL5175_probe2ensemb.csv',header = T)
结果如下
表格的第6列与第12列即为探针ID与基因的ensembID，提取这两列后将表达矩阵中的探针ID替换为ensembID，再通过org.Hs.eg.db包转换为gene_symbol即可。
网友评论