105-Bioconductor没想象的那么简单(part6)

刘小泽写于19.4.17 part1:https://www.jianshu.com/p/b030d05a80a4 part2:https://www.jianshu.com/p/c4fdab9929db part3:https://www.jianshu.com/p/f43d9d07a4af part4:https://www.jianshu.com/p/4560a29d8da2 part5:https://www.jianshu.com/p/e695b4a1dcf7 part6:https://www.jianshu.com/p/79ad7355fb8c

什么是rtracklayer?

它就是读取多种序列注释信息文件的工具,包括BED、bigwig、GTF等

比如可以在GENCODE网站下载hg38的基因组注释文件:ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_30/gencode.v30.annotation.gtf.gz

library(rtracklayer)
hg38.gencode.gtf <- import.gff(con = 'gencode.v30.annotation.gtf')
length(hg38.gencode.gtf) #2774299行数据
head(hg38.gencode.gtf) #打印前几行,还是GRange格式,注意其中gene id、transcripts id的.后内容是版本号

image-20190805151932444

还是原来的gtf文件,和linux中不同的是,它的前三列是标准GRanges格式,后面是附加的metadata

如何进行坐标转换?

有时自己分析序列差异可能是根据hg38基因组得到的,但是同时又分析了公共数据,它是基于hg19,因此为了二者可以进行比较,需要将hg38的坐标映射到hg19上。在R中需要用到chainfile文件和liftOver函数来解决

还是先加载AnnotationHub,看看其中的数据格式

library(AnnotationHub)
ah <- AnnotationHub()
table(ah$rdataclass) #其中ChainFile是需要用到的

比如想看hg38和hg19坐标转换的文件

chain <- query(ah,c("hg19","hg38","chainfile")) # 查看相关坐标文件
chain

AnnotationHub with 2 records
# snapshotDate(): 2018-10-24 
# $dataprovider: UCSC
# $species: Homo sapiens
# $rdataclass: ChainFile
# additional mcols(): taxonomyid, genome, description,
#   coordinate_1_based, maintainer, rdatadateadded, preparerclass,
#   tags, rdatapath, sourceurl, sourcetype 
# retrieve records with, e.g., 'object[["AH14108"]]' 

            title                   
  AH14108 | hg38ToHg19.over.chain.gz
  AH14150 | hg19ToHg38.over.chain.gz

可以看到hg38和hg19是可以相互转换的,只需要下载需要的文件,这里我们要将hg38转为hg19

hg38.to.hg19 <- chain[["AH14108"]] # 下载hg38转hg19坐标数据

然后开始进行坐标转换

test.hg38 <- head(hg38.gencode.gtf)
test.hg38 
test.hg19 <- liftOver(test.hg38,hg38.to.hg19) # 转换过程,先放原始数据,再放转换坐标文件
test.hg19 #得到的是一个GRangesList
new.test.hg19 <- unlist(test.hg19)

如何从GEO中下载大量的原始数据或者表达量数据?

一般一个实验会涉及十几或者几十样本,更夸张的,像单细胞数据,成白上千个单个细胞样本数据存在一个GEO ID中,虽然每个细胞的数据量不大,但是这么多如何才能下载?

以18年12月的利用smart-seq2进行的成纤维单细胞CAFs分析文章为例: Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing

它的数据存在了 GSE111229,其中有768个sample

library(GEOquery)
test.geo <- getGEO("GSE111229")

下载的这个是一个列表,用length(test.geo) 看一下就知道,结果是1,现在取出真正的数据,也就是一个ExpressionSet

> test.data <- test.geo[[1]]

ExpressionSet (storageMode: lockedEnvironment)
assayData: 0 features, 768 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: GSM3025845 GSM3025846 ... GSM3026612 (768 total)
  varLabels: title geo_accession ... tissue:ch1 (47 total)
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
  pubMedIds: 30514914 
Annotation: GPL13112 

然后expressionSet有很多参数,像上面的phenoData就可以用pData()来获取,可以用sampleNames获得样本名

pd <- pData(test.data)
sampleNames(test.data) 

pd是一个数据框,其中有一列taxid是物种分类号:10090,表示物种是小鼠;另外还有data_processing的信息,包括比对、使用RPKM定量等;再往后,relation.1中记录了每个样本的SRA位置,可以下载;有的GSE数据还给出了每个样本的表达量数据(一般在supplement data中),那么就可以$直接提取那一列,as.data.frame()保存为数据框,write.csv()输出为文本文件,最后用循环下载wget -i data.txt

Yunze Liu
Yunze Liu
Bioinformatics Sharer

Co-founder of Bioinfoplanet(生信星球)

Next
Previous

Related