差异表达edgeR，limma（上）
2019-07-03 本文已影响0人
nnlrl
第一部分 背景介绍
为了方便起见，直接选用airway的数据作为训练数据
library(airway)
library(edgeR)
library(limma)
airway
#class: RangedSummarizedExperiment
#dim: 64102 8
#metadata(1): ''
#assays(1): counts
#rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
#rowData names(0):
#colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
#colData names(9): SampleName cell ... Sample BioSample
##airway是一个`RangedSummarizedExperiment`对象，我们可以直接用`colData()`提取它的样本信息
colData(airway)
#DataFrame with 8 rows and 9 columns
# SampleName cell dex albut Run avgLength Experiment Sample BioSample
# <factor> <factor> <factor> <factor> <factor> <integer> <factor> <factor> <factor>
#SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126 SRX384345 SRS508568 SAMN02422669
#SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126 SRX384346 SRS508567 SAMN02422675
#SRR1039512 GSM1275866 N052611 untrt untrt SRR1039512 126 SRX384349 SRS508571 SAMN02422678
#SRR1039513 GSM1275867 N052611 trt untrt SRR1039513 87 SRX384350 SRS508572 SAMN02422670
#SRR1039516 GSM1275870 N080611 untrt untrt SRR1039516 120 SRX384353 SRS508575 SAMN02422682
#SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126 SRX384354 SRS508576 SAMN02422673
#SRR1039520 GSM1275874 N061011 untrt untrt SRR1039520 101 SRX384357 SRS508579 SAMN02422683
#SRR1039521 GSM1275875 N061011 trt untrt SRR1039521 98 SRX384358 SRS508580 SAMN02422677
##可以看到，第三列就是我们想要的表型信息，于是我们提取这个表型信息
pheno <- colData(airway)[,3]
##之后用assay（）提取表达矩阵
airway <- assay(airway)
head(airway)
# SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 SRR1039521
#ENSG00000000003 679 448 873 408 1138 1047 770 572
#ENSG00000000005 0 0 0 0 0 0 0 0
#ENSG00000000419 467 515 621 365 587 799 417 508
#ENSG00000000457 260 211 263 164 245 331 233 229
#ENSG00000000460 60 55 40 35 78 63 76 60
#ENSG00000000938 0 0 2 0 1 0 0 0
基因注释
##上一步我们就得到了它的表达矩阵以及表型信息，之后将对这些基因进行注释，可以看到，表达矩阵的基因都是ensembl形式的，我们要找到它的`external_gene_name`,`entrezgene`,`chromosome_name`,`gene_biotype`等注释信息已进行后续操作
library(biomaRt)
gene.names <- rownames(airway)
ensembl <- useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset="hsapiens_gene_ensembl",
host="aug2017.archive.ensembl.org") # Ensembl version 90.
ensemblGenes <- getBM(attributes=c('ensembl_gene_id', 'chromosome_name', 'gene_biotype',
'external_gene_name', 'entrezgene'), filters='ensembl_gene_id',
values=gene.names, mart=ensembl)
head(ensemblGenes)
#ensembl_gene_id chromosome_name gene_biotype external_gene_name entrezgene
#1 ENSG00000000003 X protein_coding TSPAN6 7105
#2 ENSG00000000005 X protein_coding TNMD 64102
#3 ENSG00000000419 20 protein_coding DPM1 8813
#4 ENSG00000000457 1 protein_coding SCYL3 57147
#5 ENSG00000000460 1 protein_coding C1orf112 55732
#6 ENSG00000000938 1 protein_coding FGR 2268
features <- ensemblGenes[match(gene.names, ensemblGenes$ensembl_gene_id),]
features$ensembl_gene_id <- gene.names
features$entrezgene <- gsub(" ", "", as.character(features$entrezgene))
row.names(features) <- gene.names
head(features)
# ensembl_gene_id chromosome_name gene_biotype external_gene_name entrezgene
#ENSG00000000003 ENSG00000000003 X protein_coding TSPAN6 7105
#ENSG00000000005 ENSG00000000005 X protein_coding TNMD 64102
#ENSG00000000419 ENSG00000000419 20 protein_coding DPM1 8813
#ENSG00000000457 ENSG00000000457 1 protein_coding SCYL3 57147
#ENSG00000000460 ENSG00000000460 1 protein_coding C1orf112 55732
#ENSG00000000938 ENSG00000000938 1 protein_coding FGR 2268
第二部分 得到DGEList，并对其进行操作
y = DGEList(airway)
y$samples$group <- pheno
y$features = features
##这样y中就包含了样本信息以及基因的信息
library(scater)
new.row.names <- uniquifyFeatureNames(
features$ensembl_gene_id,
features$external_gene_name
)
rownames(y) <- new.row.names
head(rownames(y), 20)
#[1] "TSPAN6" "TNMD" "DPM1" "SCYL3" "C1orf112" "FGR" "CFH" "FUCA2" "GCLC" "NFYA" "STPG1"
#[12] "NIPAL3" "LAS1L" "ENPP4" "SEMA3F" "CFTR" "ANKIB1" "CYP51A1" "KRIT1" "RAD52"
使用symbol_id为y命名，并对重复的symbol_id与ensembl_id结合以防止重复id
uniquifyFeatureNames(
ID=paste0("ENSG0000000", 1:5),
names=c("A", NA, "B", "C", "A")
)
# "A_ENSG00000001" "ENSG00000002" "B" "C" "A_ENSG00000005"
第三部分 数据预处理
原始数据的转换
进行后续处理一般都不会用raw counts的，因为存在测序深度、文库大小的差别，这样的结果是不准确的
一般的做法是：利用标准化算法，如CPM（counts per million), log-CPM (log2-counts per million), RPKM (reads per kilobase of transcript per million), FPKM(fragments per kilobase oftranscript per million)等去除文库大小、深度的影响。和RPKM、FPKM不同的是，CPM和log-CPM不需要考虑feature length的差异，也就是说基因长度在统计时被当成常数，只考虑不同处理下的不同，而不会受长度的影响
cpm使用cpm()函数；RPKM使用rpkm函数，都属于edegR
cpm <- cpm(y)
lcpm <- cpm(y, log=TRUE)
预筛选基因
所有数据集都会存在一些检测不到的基因或者表达很少的基因，这些基因对于后续分析来说不仅对结果没有影响，还会增加下游分析的计算量，所以我们需要设计一定的阈值将这些基因去除
首先看一下表达量为0的基因
table(rowSums(y$counts==0)==6)
#FALSE TRUE
#61682 2420
##在所有样本中表达量都为0的基因共有2420个
过滤基因：标准就是lcpm在所有样本表达量的平均值大于0,也就是cpm大于1
keep <- rowMeans(lcpm>0)>0
y <- y[keep ,, keep.lib.sizes=FALSE]
dim(y)
#[1] 15806 8
数据标准化
对于edgeR来说，进行差异表达时为什么需要counts矩阵呢，因为edgeR有自己的标准化步骤，edgeR使用TMM（trimmed mean of M-values）算法，利用edgeR中函数calNormFactors()
标准化用到的normalisation factors 就在DEGList中的x$samples$norm.factors
y <- calcNormFactors(y, method = "TMM")
y$samples$norm.factors
#[1] 1.0550718 1.0217435 0.9897916 0.9498495 1.0307290 0.9785444 1.0258677 0.9535888
如何检查是否标准化，可以做一个箱线图显示
#模拟一个数据
x2 <- y
x2$samples$norm.factors <- 1
x2$counts[,1] <- ceiling(x2$counts[,1]*0.05) #第一个样本count缩小到原来5%
x2$counts[,2] <- x2$counts[,2]*5 # 第二个样本扩大到原来5倍
#画图
par(mfrow=c(1,2))
lcpm <- cpm(x2, log=TRUE)
boxplot(lcpm, las=2, col=col, main="")
title(main="A. Example: Unnormalised data",ylab="Log-cpm")
x2 <- calcNormFactors(x2)
x2$samples$norm.factors
#[1] 0.06578336 6.03546029 1.16646293 1.12351485 1.21732873 1.15208812 1.21158739 1.13103354
lcpm <- cpm(x2, log=TRUE)
boxplot(lcpm, las=2, col=col, main="")
title(main="B. Example: Normalised data",ylab="Log-cpm")
箱线图查看标准化
