7 Revisión
- ¿Debemos explorar las relaciones entre nuestras variables con información de nuestras muestras previo a hacer un análisis de expresión diferencial?
- ¿Por qué usamos el paquete
edgeR? - ¿Por qué es importante el argumento
sort.byentopTable()? - ¿Por qué es importante el argumento
coefentopTable()?
Usemos los datos de http://research.libd.org/SPEAQeasy-example/bootcamp_intro
speaqeasy_data <- file.path(tempdir(), "rse_speaqeasy.RData")
download.file("https://github.com/LieberInstitute/SPEAQeasy-example/blob/master/rse_speaqeasy.RData?raw=true", speaqeasy_data, mode = "wb")
library("SummarizedExperiment")
load(speaqeasy_data, verbose = TRUE)## Loading objects:
## rse_gene## class: RangedSummarizedExperiment
## dim: 60609 40
## metadata(0):
## assays(1): counts
## rownames(60609): ENSG00000223972.5 ENSG00000227232.5 ... ENSG00000210195.2 ENSG00000210196.2
## rowData names(10): Length gencodeID ... NumTx gencodeTx
## colnames(40): R13896_H7JKMBBXX R13903_HCTYLBBXX ... R15120_HFY2MBBXX R15134_HFFGHBBXX
## colData names(67): SAMPLE_ID FQCbasicStats ... AgeDeath BrNum- ¿Cuantos genes y muestras tenemos en estos datos?
7.1 Ejercicio en equipo
- ¿Hay diferencias en
totalAssignedGeneomitoRateentre los grupos de diagnosis (PrimaryDx)? - Grafica la expresión de SNAP25 para cada grupo de diagnosis.
- Sugiere un modelo estadistico que podríamos usar en una análisis de expresión diferencial. Verifica que si sea un modelo full rank. ¿Cúal sería el o los coeficientes de interés?
7.2 Respuestas
##
## Bipolar Control Other
## 20 20 0## Eliminemos el diagnosis "Other" porque no tiene información
rse_gene$PrimaryDx <- droplevels(rse_gene$PrimaryDx)
table(rse_gene$PrimaryDx)##
## Bipolar Control
## 20 20## Exploremos numéricamente diferencias entre grupos de diagnosis para
## varias variables
with(colData(rse_gene), tapply(totalAssignedGene, PrimaryDx, summary))## $Bipolar
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.09529 0.24454 0.29773 0.28867 0.34662 0.38358
##
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1384 0.2575 0.2895 0.2984 0.3537 0.4073## $Bipolar
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.01297 0.01555 0.01735 0.01818 0.01980 0.03019
##
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.01194 0.01564 0.01911 0.02029 0.02439 0.03037## Podemos hacer lo mismo para otras variables
with(colData(rse_gene), tapply(mitoRate, BrainRegion, summary))## $Amygdala
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.01302 0.01706 0.02029 0.02097 0.02437 0.03037
##
## $sACC
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.01194 0.01449 0.01686 0.01749 0.01836 0.03019## Podemos resolver la primeras preguntas con iSEE
if (interactive()) iSEE::iSEE(rse_gene)
## O hacer graficas nosotros mismos. Aquí les muestro una posible respuesta
## con ggplot2
library("ggplot2")
ggplot(
as.data.frame(colData(rse_gene)),
aes(y = totalAssignedGene, group = PrimaryDx, x = PrimaryDx)
) +
geom_boxplot() +
theme_bw(base_size = 20) +
xlab("Diagnosis")
ggplot(
as.data.frame(colData(rse_gene)),
aes(y = totalAssignedGene, group = paste0(PrimaryDx, "_", BrainRegion), x = paste0(PrimaryDx, "_", BrainRegion))
) +
geom_boxplot() +
theme_bw(base_size = 20) +
xlab("Diagnosis")
ggplot(
as.data.frame(colData(rse_gene)),
aes(y = mitoRate, group = PrimaryDx, x = PrimaryDx)
) +
geom_boxplot() +
theme_bw(base_size = 20) +
xlab("Diagnosis")
ggplot(
as.data.frame(colData(rse_gene)),
aes(y = mitoRate, group = paste0(PrimaryDx, "_", BrainRegion), x = paste0(PrimaryDx, "_", BrainRegion))
) +
geom_boxplot() +
theme_bw(base_size = 20) +
xlab("Diagnosis")
## Otras variables
ggplot(
as.data.frame(colData(rse_gene)),
aes(y = mitoRate, group = BrainRegion, x = BrainRegion)
) +
geom_boxplot() +
theme_bw(base_size = 20) +
xlab("Brain Region")
## GRanges object with 60609 ranges and 10 metadata columns:
## seqnames ranges strand | Length gencodeID ensemblID gene_type
## <Rle> <IRanges> <Rle> | <integer> <character> <character> <character>
## ENSG00000223972.5 chr1 11869-13670 + | 1735 ENSG00000223972.5 ENSG00000223972 transcribed_unproces..
## ENSG00000227232.5 chr1 14404-29570 - | 1351 ENSG00000227232.5 ENSG00000227232 unprocessed_pseudogene
## ENSG00000278267.1 chr1 17369-17436 - | 68 ENSG00000278267.1 ENSG00000278267 miRNA
## ENSG00000243485.5 chr1 29554-31097 + | 1021 ENSG00000243485.5 ENSG00000243485 lncRNA
## ENSG00000284332.1 chr1 30366-30503 + | 138 ENSG00000284332.1 ENSG00000284332 miRNA
## ... ... ... ... . ... ... ... ...
## ENSG00000198695.2 chrM 14149-14673 - | 525 ENSG00000198695.2 ENSG00000198695 protein_coding
## ENSG00000210194.1 chrM 14674-14742 - | 69 ENSG00000210194.1 ENSG00000210194 Mt_tRNA
## ENSG00000198727.2 chrM 14747-15887 + | 1141 ENSG00000198727.2 ENSG00000198727 protein_coding
## ENSG00000210195.2 chrM 15888-15953 + | 66 ENSG00000210195.2 ENSG00000210195 Mt_tRNA
## ENSG00000210196.2 chrM 15956-16023 - | 68 ENSG00000210196.2 ENSG00000210196 Mt_tRNA
## Symbol EntrezID Class meanExprs NumTx gencodeTx
## <character> <integer> <character> <numeric> <integer> <character>
## ENSG00000223972.5 DDX11L1 <NA> InGen 0.00347263 2 ENST00000456328.2;EN..
## ENSG00000227232.5 WASH7P <NA> InGen 1.57167877 1 ENST00000488147.1
## ENSG00000278267.1 MIR6859-1 <NA> InGen 7.74437941 1 ENST00000619216.1
## ENSG00000243485.5 MIR1302-2HG <NA> InGen 0.01403484 2 ENST00000473358.1;EN..
## ENSG00000284332.1 MIR1302-2 <NA> InGen 0.00000000 1 ENST00000607096.1
## ... ... ... ... ... ... ...
## ENSG00000198695.2 MT-ND6 4541 InGen 182.92863 1 ENST00000361681.2
## ENSG00000210194.1 MT-TE <NA> InGen 9.80684 1 ENST00000387459.1
## ENSG00000198727.2 MT-CYB 4519 InGen 826.87424 1 ENST00000361789.2
## ENSG00000210195.2 MT-TT <NA> InGen 276.37225 1 ENST00000387460.2
## ENSG00000210196.2 MT-TP <NA> InGen 366.58138 1 ENST00000387461.2
## -------
## seqinfo: 25 sequences from an unspecified genome; no seqlengths## En este objeto los nombres de los genes vienen en la variable "Symbol"
i <- which(rowRanges(rse_gene)$Symbol == "SNAP25")
i## [1] 54098## Para graficar con ggplot2, hagamos un pequeño data.frame
df <- data.frame(
expression = assay(rse_gene)[i, ],
Dx = rse_gene$PrimaryDx
)
## Ya teniendo el pequeño data.frame, podemos hacer la gráfica
ggplot(df, aes(y = log2(expression + 0.5), group = Dx, x = Dx)) +
geom_boxplot() +
theme_bw(base_size = 20) +
xlab("Diagnosis") +
ylab("SNAP25: log2(x + 0.5)")
## https://bioconductor.org/packages/release/bioc/vignettes/scater/inst/doc/overview.html#3_Visualizing_expression_values
scater::plotExpression(
as(rse_gene, "SingleCellExperiment"),
features = rownames(rse_gene)[i],
x = "PrimaryDx",
exprs_values = "counts",
colour_by = "BrainRegion",
xlab = "Diagnosis"
)
if (requireNamespace("plotly", quietly = TRUE)) {
## Lo pueden instalar con
# install.packages("plotly")
## Guardemos el resultado de plotExpression()
p <- scater::plotExpression(
as(rse_gene, "SingleCellExperiment"),
features = rownames(rse_gene)[i],
x = "PrimaryDx",
exprs_values = "counts",
colour_by = "BrainRegion",
xlab = "Diagnosis"
)
## scater::plotExpression() regresa un objeto de clase ggplot
class(p)
## así que podemos usar plotly para crear una versión
## interactiva
plotly::ggplotly(p)
}## [1] "SAMPLE_ID" "FQCbasicStats" "perBaseQual"
## [4] "perTileQual" "perSeqQual" "perBaseContent"
## [7] "GCcontent" "Ncontent" "SeqLengthDist"
## [10] "SeqDuplication" "OverrepSeqs" "AdapterContent"
## [13] "KmerContent" "SeqLength_R1" "percentGC_R1"
## [16] "phred20-21_R1" "phred48-49_R1" "phred76-77_R1"
## [19] "phred100-101_R1" "phredGT30_R1" "phredGT35_R1"
## [22] "Adapter50-51_R1" "Adapter70-71_R1" "Adapter90_R1"
## [25] "SeqLength_R2" "percentGC_R2" "phred20-21_R2"
## [28] "phred48-49_R2" "phred76-77_R2" "phred100-101_R2"
## [31] "phredGT30_R2" "phredGT35_R2" "Adapter50-51_R2"
## [34] "Adapter70-71_R2" "Adapter90_R2" "bamFile"
## [37] "trimmed" "numReads" "numMapped"
## [40] "numUnmapped" "overallMapRate" "concordMapRate"
## [43] "totalMapped" "mitoMapped" "mitoRate"
## [46] "totalAssignedGene" "gene_Assigned" "gene_Unassigned_Unmapped"
## [49] "gene_Unassigned_Read_Type" "gene_Unassigned_Singleton" "gene_Unassigned_MappingQuality"
## [52] "gene_Unassigned_Chimera" "gene_Unassigned_FragmentLength" "gene_Unassigned_Duplicate"
## [55] "gene_Unassigned_MultiMapping" "gene_Unassigned_Secondary" "gene_Unassigned_NonSplit"
## [58] "gene_Unassigned_NoFeatures" "gene_Unassigned_Overlapping_Length" "gene_Unassigned_Ambiguity"
## [61] "rRNA_rate" "BrainRegion" "Race"
## [64] "PrimaryDx" "Sex" "AgeDeath"
## [67] "BrNum"##
## Amygdala sACC
## 20 20##
## CAUC
## 40## Ojo! Acá es importante que hayamos usado droplevels(rse_gene$PrimaryDx)
## si no, vamos a tener un modelo que no sea _full rank_
mod <- with(
colData(rse_gene),
model.matrix(~ PrimaryDx + totalAssignedGene + mitoRate + rRNA_rate + BrainRegion + Sex + AgeDeath)
)
## Exploremos el modelo de forma interactiva
if (interactive()) {
## Tenemos que eliminar columnas que tienen NAs.
info_no_NAs <- colData(rse_gene)[, c(
"PrimaryDx", "totalAssignedGene", "rRNA_rate", "BrainRegion", "Sex",
"AgeDeath", "mitoRate", "Race"
)]
ExploreModelMatrix::ExploreModelMatrix(
info_no_NAs,
~ PrimaryDx + totalAssignedGene + mitoRate + rRNA_rate + BrainRegion + Sex + AgeDeath
)
## Veamos un modelo más sencillo sin las variables numéricas (continuas) porque
## ExploreModelMatrix nos las muestra como si fueran factors (categoricas)
## en vez de continuas
ExploreModelMatrix::ExploreModelMatrix(
info_no_NAs,
~ PrimaryDx + BrainRegion + Sex
)
## Si agregamos + Race nos da errores porque Race solo tiene 1 opción
# ExploreModelMatrix::ExploreModelMatrix(
# info_no_NAs,
# ~ PrimaryDx + BrainRegion + Sex + Race
# )
}¿Quieres más datos? Tenemos muchos en LIBD incluyendo http://eqtl.brainseq.org/phase2/.