Multimodal analysis
Data preparation
Load libraries
Load the libraries required.
library(RMINC)## Loading required package: BatchJobs## Loading required package: BBmisc## Sourcing configuration file: '/axiom2/projects/software/arch/linux-xenial-xerus/RMINC/1.5.0.0/BatchJobs/etc/BatchJobs_global_config.R'## BatchJobs configuration:
## cluster functions: Interactive
## mail.from:
## mail.to:
## mail.start: none
## mail.done: none
## mail.error: none
## default.resources:
## debug: FALSE
## raise.warnings: FALSE
## staged.queries: TRUE
## max.concurrent.jobs: Inf
## fs.timeout: NA## Sourcing configuration file: '/axiom2/projects/software/arch/linux-xenial-xerus/RMINC/1.5.0.0/RMINC/parallel/sge_BatchJobs.R'library(data.tree)
library(ggplot2)Load required data
First we point to the location of the data and atlases.
datafile<-"/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/stats/datatable_multimodal.csv"
atlasLabels<-"/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/stats/DSU_v2_mapping_of_labels_withABI.csv"
resampAtlas<-"/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/atlas/DSU_v2_on_NRXN1a_labels_resampled.mnc"
resampMask <- "/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/atlas/Dorr_2008_Steadman_2013_Ullmann_2013_on_NRXN1a_v1_mask_resampled.mnc"
cortexMask <- "/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/stats/mask_cortex.mnc"Next, read in the data.
defs <- read.csv(atlasLabels)
df <- as.data.frame(read.csv(datafile, header = TRUE, sep = ",", stringsAsFactors=FALSE), stringsAsFactors=FALSE)Get irradiation dataset brain region volumes
Now, we want to extract brain region volumes from a set of image determinant (MINC) files. The regions are defined in the atlas; the location of this atlas is in resampAtlas. The anatGetAll() function from our RMINC package will compute the volumes under these labels
vols <- anatGetAll(filename = df$log_det_scaled0.1, atlas = resampAtlas, defs=atlasLabels, method="jacobians") #Find volumes of structures from resampled atlas
volsCombined <- anatCombineStructures(vols, defs=atlasLabels, method="jacobians") #Get volumes of the combined (left and right) structuresGet MTR metric
Next, do the same thing to get the mean MTR value under each label. Note that this is not volumetric data, and so the method argument must be accordingly specified.
mtr <- anatGetAll(filename = df$mtr_orig, atlas = resampAtlas, defs=atlasLabels, method="means")
mtrCombined <- anatCombineStructures(mtr, defs=atlasLabels, method="means")
mtrCortex <- anatGetAll(filename = df$mtr_orig, atlas = cortexMask, defs=atlasLabels, method="means")## Missing Structures found in labels but not in any files: right amygdala, right Amygdalopiriform transition area, right anterior commissure: pars anterior, right anterior commissure: pars posterior, right anterior lobule (lobules 4-5), right anterior lobule white matter, right basal forebrain, right bed nucleus of stria terminalis, right Caudomedial entorhinal cortex, right cerebellar peduncle: inferior, right cerebellar peduncle: middle, right cerebellar peduncle: superior, right cerebral aqueduct, right cerebral peduncle, right Cingulate cortex: area 24a, right Cingulate cortex: area 24a', right Cingulate cortex: area 24b, right Cingulate cortex: area 24b', right Cingulate cortex: area 25, right Cingulate cortex: area 29a, right Cingulate cortex: area 29b, right Cingulate cortex: area 29c, right Cingulate cortex: area 30, right Cingulate cortex: area 32, right Cingulum , right Claustrum, right Claustrum: dorsal part, right Claustrum: ventral part, right colliculus: inferior, right colliculus: superior, right copula white matter, right copula: pyramis (lobule 8), right corpus callosum, right Cortex-amygdala transition zones, right corticospinal tract/pyramids, right crus 1 white matter, right crus 1: ansiform lobule (lobule 6), right crus 2 white matter, right crus 2: ansiform lobule (lobule 7), right cuneate nucleus, right dentate gyrus of hippocampus, right Dorsal intermediate entorhinal cortex, right Dorsal nucleus of the endopiriform, right Dorsal tenia tecta, right Dorsolateral entorhinal cortex, right Dorsolateral orbital cortex, right Ectorhinal cortex, right facial nerve (cranial nerve 7), right fasciculus retroflexus, right fastigial nucleus, right fimbria, right flocculus (FL), right flocculus white matter, right fornix, right fourth ventricle, right Frontal association cortex, right Frontal cortex: area 3, right fundus of striatum, right globus pallidus, right habenular commissure, right hippocampus, right hypothalamus, right inferior olivary complex, right Insular region: not subdivided, right Intermediate nucleus of the endopiriform claustrum, right internal capsule, right interpedunclar nucleus, right lateral olfactory tract, right Lateral orbital cortex, right Lateral parietal association cortex, right lateral septum, right lateral ventricle, right lobule 1-2 white matter, right lobule 10 white matter, right lobule 10: nodulus, right lobule 3 white matter, right lobule 3: central lobule (dorsal), right lobule 6: declive, right lobule 7: tuber (or folium), right lobule 8 white matter, right lobule 8: pyramis, right lobule 9 white matter, right lobule 9: uvula, right lobules 1-2: lingula and central lobule (ventral), right lobules 4-5 white matter, right lobules 4-5: culmen (ventral and dorsal), right lobules 6-7 white matter, right mammillary bodies, right mammilothalamic tract, right Medial entorhinal cortex, right medial lemniscus/medial longitudinal fasciculus, right Medial orbital cortex, right Medial parietal association cortex, right medial septum, right medulla, right midbrain, right nucleus accumbens, right nucleus interpositus, right olfactory bulbs, right olfactory tubercle, right optic tract, right paraflocculus (PFL), right paraflocculus white matter, right paramedian lobule, right paramedian lobule (lobule 7), right Parietal cortex: posterior area: rostral part, right periaqueductal grey, right Perirhinal cortex, right Piriform cortex, right pons, right pontine nucleus, right posterior commissure, right Posterolateral cortical amygdaloid area , right Posteromedial cortical amygdaloid area , right pre-para subiculum, right Primary auditory cortex, right Primary motor cortex, right Primary somatosensory cortex, right Primary somatosensory cortex: barrel field, right Primary somatosensory cortex: dysgranular zone, right Primary somatosensory cortex: forelimb region, right Primary somatosensory cortex: hindlimb region, right Primary somatosensory cortex: jaw region, right Primary somatosensory cortex: shoulder region, right Primary somatosensory cortex: trunk region, right Primary somatosensory cortex: upper lip region, right Primary visual cortex, right Primary visual cortex: binocular area, right Primary visual cortex: monocular area, right Rostral amygdalopiriform area, right Secondary auditory cortex: dorsal area, right Secondary auditory cortex: ventral area, right Secondary motor cortex, right Secondary somatosensory cortex, right Secondary visual cortex: lateral area, right Secondary visual cortex: mediolateral area, right Secondary visual cortex: mediomedial area, right simple lobule (lobule 6), right simple lobule white matter, right stratum granulosum of hippocampus, right stria medullaris, right stria terminalis, right striatum, right subependymale zone / rhinocele, right superior olivary complex, right Temporal association area, right thalamus, right third ventricle, right trunk of arbor vita, right trunk of crus 2 and paramedian white matter, right trunk of lobules 1-3 white matter, right trunk of lobules 6-8 white matter, right trunk of simple and crus 1 white matter, right Ventral intermediate entorhinal cortex, right Ventral nucleus of the endopiriform claustrum, right Ventral orbital cortex, right ventral tegmental decussation, right Ventral tenia tecta, left amygdala, left Amygdalopiriform transition area, left anterior commissure: pars anterior, left anterior commissure: pars posterior, left anterior lobule (lobules 4-5), left anterior lobule white matter, left basal forebrain, left bed nucleus of stria terminalis, left Caudomedial entorhinal cortex, left cerebellar peduncle: inferior, left cerebellar peduncle: middle, left cerebellar peduncle: superior, left cerebral aqueduct, left cerebral peduncle, left Cingulate cortex: area 24a, left Cingulate cortex: area 24a', left Cingulate cortex: area 24b, left Cingulate cortex: area 24b', left Cingulate cortex: area 25, left Cingulate cortex: area 29a, left Cingulate cortex: area 29b, left Cingulate cortex: area 29c, left Cingulate cortex: area 30, left Cingulate cortex: area 32, left Cingulum , left Claustrum, left Claustrum: dorsal part, left Claustrum: ventral part, left colliculus: inferior, left colliculus: superior, left copula white matter, left copula: pyramis (lobule 8), left corpus callosum, left Cortex-amygdala transition zones, left corticospinal tract/pyramids, left crus 1 white matter, left crus 1: ansiform lobule (lobule 6), left crus 2 white matter, left crus 2: ansiform lobule (lobule 7), left cuneate nucleus, left dentate gyrus of hippocampus, left dentate nucleus, left Dorsal intermediate entorhinal cortex, left Dorsal nucleus of the endopiriform, left Dorsal tenia tecta, left Dorsolateral entorhinal cortex, left Dorsolateral orbital cortex, left Ectorhinal cortex, left facial nerve (cranial nerve 7), left fasciculus retroflexus, left fastigial nucleus, left fimbria, left flocculus (FL), left flocculus white matter, left fornix, left fourth ventricle, left Frontal association cortex, left Frontal cortex: area 3, left fundus of striatum, left globus pallidus, left habenular commissure, left hippocampus, left hypothalamus, left inferior olivary complex, left Insular region: not subdivided, left Intermediate nucleus of the endopiriform claustrum, left internal capsule, left interpedunclar nucleus, left lateral olfactory tract, left Lateral orbital cortex, left Lateral parietal association cortex, left lateral septum, left lateral ventricle, left lobule 1-2 white matter, left lobule 10 white matter, left lobule 10: nodulus, left lobule 3 white matter, left lobule 3: central lobule (dorsal), left lobule 6: declive, left lobule 7: tuber (or folium), left lobule 8 white matter, left lobule 8: pyramis, left lobule 9 white matter, left lobule 9: uvula, left lobules 1-2: lingula and central lobule (ventral), left lobules 4-5 white matter, left lobules 4-5: culmen (ventral and dorsal), left lobules 6-7 white matter, left mammillary bodies, left mammilothalamic tract, left Medial entorhinal cortex, left medial lemniscus/medial longitudinal fasciculus, left Medial orbital cortex, left Medial parietal association cortex, left medial septum, left medulla, left midbrain, left nucleus accumbens, left nucleus interpositus, left olfactory bulbs, left olfactory tubercle, left optic tract, left paraflocculus (PFL), left paraflocculus white matter, left paramedian lobule, left paramedian lobule (lobule 7), left Parietal cortex: posterior area: rostral part, left periaqueductal grey, left Perirhinal cortex, left Piriform cortex, left pons, left pontine nucleus, left posterior commissure, left Posterolateral cortical amygdaloid area , left Posteromedial cortical amygdaloid area , left pre-para subiculum, left Primary auditory cortex, left Primary motor cortex, left Primary somatosensory cortex, left Primary somatosensory cortex: barrel field, left Primary somatosensory cortex: dysgranular zone, left Primary somatosensory cortex: forelimb region, left Primary somatosensory cortex: hindlimb region, left Primary somatosensory cortex: jaw region, left Primary somatosensory cortex: shoulder region, left Primary somatosensory cortex: trunk region, left Primary somatosensory cortex: upper lip region, left Primary visual cortex, left Primary visual cortex: binocular area, left Primary visual cortex: monocular area, left Rostral amygdalopiriform area, left Secondary auditory cortex: dorsal area, left Secondary auditory cortex: ventral area, left Secondary motor cortex, left Secondary somatosensory cortex, left Secondary visual cortex: lateral area, left Secondary visual cortex: mediolateral area, left Secondary visual cortex: mediomedial area, left simple lobule (lobule 6), left simple lobule white matter, left stratum granulosum of hippocampus, left stria medullaris, left stria terminalis, left striatum, left subependymale zone / rhinocele, left superior olivary complex, left Temporal association area, left thalamus, left third ventricle, left trunk of arbor vita, left trunk of crus 2 and paramedian white matter, left trunk of lobules 1-3 white matter, left trunk of lobules 6-8 white matter, left trunk of simple and crus 1 white matter, left Ventral intermediate entorhinal cortex, left Ventral nucleus of the endopiriform claustrum, left Ventral orbital cortex, left ventral tegmental decussation, left Ventral tenia tectadf$mtrCortex <- mtrCortex[,1]Get DTI metrics
Repeat for the mean diffusivity (MD)… ### Mean diffusivity
md <- anatGetAll(filename = df$MD_orig, atlas = resampAtlas,defs=atlasLabels,method="means")
mdCombined<-anatCombineStructures(md,defs=atlasLabels,method="means")…and the fractional anistropy (FA). ### Fractional anisotropy
fa <- anatGetAll(filename = df$FA_orig, atlas = resampAtlas,defs=atlasLabels,method="means")
faCombined<-anatCombineStructures(fa,defs=atlasLabels,method="means")Analysis
MTR
Here, we will set up a linear model to model the effects of dose on MTR, for each voxel in the brain. The extra term mtrCortex is to normalize to the mean value of the MTR in the cortex.
mtr_lm <- mincLm(mtr_blur_0.2 ~ Dose + mtrCortex, df, mask=resampMask)## Method: lm
## Number of volumes: 15
## Volume sizes: 152 320 225
## N: 15 P: 3
## In slice
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151
## DoneLet’s look at what the distribution of voxel values in the average MTR image looks like…
anatFile <- "/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/transformed_mnc_files/mt/mtr_average.mnc"
anatVolume <- mincArray(mincGetVolume(anatFile))
hist(anatVolume)
…and plot this as brain slices.
mincPlotSliceSeries(anatVolume, mincArray(mtr_lm, "tvalue-Dose"), anatLow=0, anatHigh=0.2, low=2, high=15, symmetric=T, begin=50, end=-40, legend="t-statistic", plottitle = "Dose (7 vs 0 only)")
Next, compute thresholds using the false discovery rate method described by Benjamini and Hochberg (1995).
qvs <- mincFDR(mtr_lm, mask=resampMask)##
## Computing FDR threshold for all columns
## Computing threshold for F-statistic
## Computing threshold for tvalue-(Intercept)
## Computing threshold for tvalue-Dose
## Computing threshold for tvalue-mtrCortexLet’s see the thresholds:
qvs## Multidimensional MINC volume
## Columns: F-statistic tvalue-(Intercept) tvalue-Dose tvalue-mtrCortex
## [1] "/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/transformed_mnc_files/mt/10mar13_mtr.16.jan11_dc-resampled_fwhm_0.2_blur.mnc"
## Degrees of Freedom: c(2, 12) 12 12 12
## FDR Thresholds:
## F-statistic tvalue-(Intercept) tvalue-Dose tvalue-mtrCortex
## 0.01 7.445056 NA 5.289509 3.163634
## 0.05 4.059829 NA 3.612063 2.241895
## 0.1 2.910550 NA 2.879874 1.832613
## 0.15 2.305699 NA 2.472932 1.581242
## 0.2 1.902976 NA 2.180855 1.394536And let’s replot this with the colours set to reflect our FDR method determined thresholds.
mincPlotSliceSeries(anatVolume, mincArray(mtr_lm, "tvalue-Dose"), anatLow=0, anatHigh=0.2, low=3.6, high=15, symmetric=T, begin=50, end=-40, legend="t-statistic", plottitle = "Dose (7 vs 0 only)")
FA
Repeat all the above for FA.
fa_lm <- mincLm(FA_blur_0.2 ~ Dose, df, mask=resampMask)## Method: lm
## Number of volumes: 15
## Volume sizes: 152 320 225
## N: 15 P: 2
## In slice
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151
## DoneanatFile <- "/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/transformed_mnc_files/dti/FA_average.mnc"
anatVolume <- mincArray(mincGetVolume(anatFile))
hist(anatVolume)
mincPlotSliceSeries(anatVolume, mincArray(fa_lm , "tvalue-Dose"), anatLow=0, anatHigh=0.8, low=2, high=15, symmetric=T, begin=50, end=-40, legend="t-statistic", plottitle = "Dose (7 vs 0 only)")
qvs <- mincFDR(fa_lm, mask=resampMask)##
## Computing FDR threshold for all columns
## Computing threshold for F-statistic
## Computing threshold for tvalue-(Intercept)
## Computing threshold for tvalue-Doseqvs## Multidimensional MINC volume
## Columns: F-statistic tvalue-(Intercept) tvalue-Dose
## [1] "/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/transformed_mnc_files/dti/21may13.2_FA.nov07_dc-resampled_fwhm_0.2_blur.mnc"
## Degrees of Freedom: c(1, 13) 13 13
## FDR Thresholds:
## F-statistic tvalue-(Intercept) tvalue-Dose
## 0.01 NA 3.260863 NA
## 0.05 NA 3.260863 NA
## 0.1 NA 3.260863 NA
## 0.15 NA 3.260863 NA
## 0.2 84.23801 3.260863 9.178127mincPlotSliceSeries(anatVolume, mincArray(fa_lm, "tvalue-Dose"), anatLow=0, anatHigh=0.8, low=9.18, high=15, symmetric=T, begin=50, end=-40, legend="t-statistic", plottitle = "Dose (7 vs 0 only)")
Compare FA and MTR
Do FA and MTR show similar changes after irradiation? Let’s plot the t-statistic value for FA changes against the t-stat value for MTR changes. To do so, first we will determine which voxels are within the brain using the mask, and then choose only those voxels.
maskdata <- mincGetVolume(resampMask)
fa_voxels <- fa_lm[,"tvalue-Dose"][maskdata > 0.5]
mtr_voxels <- mtr_lm[,"tvalue-Dose"][maskdata > 0.5]Now for the plotting, with the ggplot2 library.
mtr_fa_voxels <- data.frame(mtr=mtr_voxels, fa=fa_voxels)
plt <- ggplot(data = mtr_fa_voxels, aes(x=mtr, y=fa))
plt <- plt + geom_point(size=0.1, alpha=0.005)
print(plt)
Let’s visualize the spatial MTR and FA changes by gridding the above plot, and labeling each voxel with a value based on which part of the grid the datapoint lies within.
mtrp5 <- mtr_fa_voxels$mtr >= 5
mtr0 <- (mtr_fa_voxels$mtr < 5 & mtr_fa_voxels$mtr >= -5)
mtrn5 <- mtr_fa_voxels$mtr < -5
fap3 <- mtr_fa_voxels$fa >= 3
fa0 <- (mtr_fa_voxels$fa < 3 & mtr_fa_voxels$fa >= -3)
fan3 <- mtr_fa_voxels$fa < -3
mtr_fa_labeled_data <- maskdata
mtr_fa_labeled_data[maskdata > 0.5][which(mtrp5 & fap3)] <- 1
mtr_fa_labeled_data[maskdata > 0.5][which(mtrp5 & fa0)] <- 2
mtr_fa_labeled_data[maskdata > 0.5][which(mtrp5 & fan3)] <- 3
mtr_fa_labeled_data[maskdata > 0.5][which(mtr0 & fap3)] <- 4
mtr_fa_labeled_data[maskdata > 0.5][which(mtr0 & fa0)] <- 0
mtr_fa_labeled_data[maskdata > 0.5][which(mtr0 & fan3)] <- 6
mtr_fa_labeled_data[maskdata > 0.5][which(mtrn5 & fap3)] <- 7
mtr_fa_labeled_data[maskdata > 0.5][which(mtrn5 & fa0)] <- 8
mtr_fa_labeled_data[maskdata > 0.5][which(mtrn5 & fan3)] <- 9
mincPlotSliceSeries(anatomy = anatVolume, mincArray(mtr_fa_labeled_data), anatLow=0, anatHigh=0.8, low=0, high=9, symmetric=F, begin=50, end=-40, legend="Clusters", plottitle = "MTR and FA ...", col = rainbow(10))
# Write the data
mincWriteVolume(mtr_fa_labeled_data, output.filename = "mtr_fa_clusters.mnc")## Warning: the outputfile already exists, continue? (y/n)
## Range: 9.000000 0.000000Gene enrichment analysis with external dataset
Now, let’s see which genes share a similar expression profile to the MTR values. First, load in the data which already has been collected for you.
load("/hpf/largeprojects/MICe/yyee/summerschool/multimodal_analysis/demo/vols_gene_and_connectivity.RData")Darren Fernandes has provided code to work with raw Allen Institute data: https://github.com/DJFernandes/ABIgeneRMINC .
The expression dataset needs a bit of preprocessing. For example, the raw values reported in the expression datasets come from different experiments, and so are not comparable. To get around this, we need to normalize the data. We’ll do a crude normalization by dividing by the total expression (row sum).
gs <- dat$expression_sagittal
gs <- gs/rowSums(gs, na.rm = TRUE)Next, only choose the structures in which the data was collected (left hemisphere)
strucs_right <- colnames(gs)[which(startsWith(colnames(gs), "right"))]
strucs_left <- colnames(gs)[which(startsWith(colnames(gs), "left"))]
strucs <- strucs_left
vl <- gs[,which(colnames(gs) %in% strucs)]Let’s visualize the distribution of relative expression levels under each structure. In particular, notice the outliers on the left end of the distribution.
hist(colSums(vl, na.rm=TRUE), 20)
Digging a bit deeper, we can find out which structures these are.
names(which(colSums(vl, na.rm = TRUE) < 30))## [1] "left Primary auditory cortex"
## [2] "left Secondary auditory cortex: ventral area"
## [3] "left Dorsolateral entorhinal cortex"
## [4] "left Ectorhinal cortex"
## [5] "left Perirhinal cortex"
## [6] "left Temporal association area"You might notice these are structures located at the most lateral part of the brain. This makes sense; the data was acquired by slicing sagitally, which might cause issues for these structures. Let’s remove them, along with CSF structures.
strucs <- setdiff(strucs, names(which(colSums(vl, na.rm = TRUE) < 30)))
strucs <- setdiff(strucs, paste("left", as.character(defs$Structure[which(defs$tissue.type =="CSF")])))
vl <- vl[,strucs]The Allen Institute also provides a list of non-expressing genes (supplementary table in their 2007 paper by Lein et. al.), so let’s also remove these. To do so, we will first read in this list. Then we’ll rename the rows of the expression dataframe to strip out the expression ID and just keep the gene name part. Lastly, we’ll subset the data to only have the rows with gene names not in the non-expressing gene list.
nonexpressing_list_of_genes <- "/hpf/largeprojects/MICe/yyee/projects/greedy_gene_enrichment/nonexpressing_genes.txt"
gene_ids_raw <- rownames(vl)
gene_ids <- gsub("\\*", "", gene_ids_raw)
gene_names <- sapply(strsplit(gene_ids, "_sid"), "[[", 1)
nonexpressors <- as.character(read.table(nonexpressing_list_of_genes)$V1)
expressors <- setdiff(gene_names, nonexpressors)
vl <- vl[which(gene_names %in% expressors),]
rownames(vl) <- gsub("\\*", "", rownames(vl))Finally, remove any genes that are not in the list provided but still have no expression.
vl[is.na(vl)] <- 0
vl <- vl[which(rowSums(vl) > 0.2),]
colnames(vl) <- gsub("left ", "", colnames(vl))Now we have a g x p table which represents the relative spatial expression pattern for g genes and p regions. We also have a structural imaging dataset (MTR) that is an n x p table, where n is the number of mice and p is the number of regions. Let’s just pick the mean MTR values for a single mouse, to give us a length p vector. We can ask the question, which genes show a similar expression profile as this MTR vector? To determine this, we will compute univariate correlations between the MTR vector and the spatial expression profile of each gene. Note that there are other methods that might be more appropriate for picking genes that associate with a structural profile such as sparse regression (LASSO, elastic nets), but let’s stick with the good old Pearson correlation for now.
mtr_1mouse <- mtrCombined[1,colnames(vl)]
cordf <- data.frame(gene_id=rownames(vl), gene_mtr_cor=numeric(dim(vl)[1]))
for (i in 1:dim(cordf)[1]) {
r <- cor(mtr_1mouse, as.numeric(vl[i,]))
cordf$gene_mtr_cor[i] <- r
#cat(".")
}
cat("Done.\n")## Done.Write out the ranked gene list to a file. If you are confused with the long expression given to the write.table function, run the nested parts of it, bit by bit, to see what’s happening.
write.table(sapply(strsplit(as.character(cordf[order(cordf$gene_mtr_cor, decreasing = TRUE),"gene_id"]), "_sid"), "[[", 1), "/tmp/MTR_1mouse_genes.txt", row.names = FALSE, col.names = F, quote=F)You can input this gene list into GO enrichment software. Here is the result for these genes: http://cbl-gorilla.cs.technion.ac.il/GOrilla/33gwj7b4/GOResults.html
Connectivity analysis
Let’s get an idea of what the connectivity data looks like.
tracer_image <- "/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/stats/data/projection/180436360/projection_density.mnc"
tracer_vol <- mincArray(mincGetVolume(tracer_image))
anatFile <- "/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/transformed_mnc_files/mt/mtr_average.mnc"
anatVolume <- mincArray(mincGetVolume(anatFile))
mincPlotSliceSeries(anatVolume, tracer_vol, anatLow=0, anatHigh=0.2, low=0.1, high=0.5, symmetric=F, begin=50, end=-40, legend="tracer density", plottitle = "", col=colorRampPalette(c("black", "green"))(100))
Get injection experiment data.
injection_experiments <- read.csv("/hpf/largeprojects/MICe/yyee/analysis/structural_covariance/voxelwise/injection_experiments.csv")[,c(1:19,22:23)]Let us try and see if the irradiation-induced changes in MD are between regions connected by fiber tracts. First, we’ll run a linear model for MD and plot the MD changes. For simplicity, let’s only consider irradiation includes decreases in MD.
md_lm <- mincLm(MD_blur_0.2 ~ Dose, df, mask=resampMask)## Method: lm
## Number of volumes: 15
## Volume sizes: 152 320 225
## N: 15 P: 2
## In slice
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151
## DoneanatFile <- "/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/transformed_mnc_files/mt/mtr_average.mnc"
anatVolume <- mincArray(mincGetVolume(anatFile))
mincPlotSliceSeries(anatVolume, mincArray(md_lm, "tvalue-Dose"), anatLow=0, anatHigh=0.2, low=2, high=15, symmetric=T, begin=50, end=-40, legend="t-statistic", plottitle = "Dose (7 vs 0 only)")
We’ve binarized and labeled the regions of negative changes in MD for you. Load that in and visualize it.
neg_labels <- "/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/stats/md_strucatlas.mnc"
neg_labels_array <- mincArray(mincGetVolume(neg_labels))
max_label_value <- max(neg_labels_array)
mincPlotSliceSeries(anatVolume, neg_labels_array, anatLow=0, anatHigh=0.2, low=1, high=max_label_value, symmetric=F, begin=50, end=-40, legend="Label", plottitle = "Dose (7 vs 0 only)", col=rainbow(151))
Let’s get the mean tracer projection value under each of these labels.
projection_directory <- "/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/stats/data/projection"
tracer_ids <- list.files(projection_directory)
tracer_files <- paste(projection_directory, tracer_ids, "projection_density.mnc", sep="/")
defs_md_neg_blobs <- "/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/stats/neg_blob_defs.csv"
#tracer_vols <- anatGetAll(tracer_files, atlas=neg_labels, defs=defs_md_neg_blobs, method="means")
#save("tracer_vols", file = "/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/stats/tracer_volumes_under_md_blobs.RData")
load("/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/stats/tracer_volumes_under_md_blobs.RData")Now, let’s count, for each tracer (row), the number of structures that have a tracer mean value of greater than 0.1.
num_elements_greater_than_0.1 <- function(x) {
return(length(which(x > 0.1)))
}
num_strucs <- apply(tracer_vols, 1, FUN = "num_elements_greater_than_0.1")
tracer_overlaps <- data.frame(tracer_id=tracer_ids, num_overlapping_strucs=num_strucs)
tracer_overlaps[which.max(tracer_overlaps$num_overlapping_strucs),]## tracer_id num_overlapping_strucs
## 339 158257355 3injection_experiments[injection_experiments$id==tracer_overlaps$tracer_id[which.max(tracer_overlaps$num_overlapping_strucs)],]## X id name structure.id structure.abbrev
## 7 10 158257355 378-1626 596 NDB
## structure.name structure.color strain transgenic.line gender
## 7 Diagonal band nucleus 96A7D3 C57BL/6J NA M
## injection.coordinates injection.volume
## 7 [4440, 6130, 5460] 1.41606
## injection.structures
## 7 [{"id"=>48, "name"=>"Anterior cingulate area, ventral part", "abbreviation"=>"ACAv", "color"=>"40A666"}, {"id"=>56, "name"=>"Nucleus accumbens", "abbreviation"=>"ACB", "color"=>"80CDF8"}, {"id"=>72, "name"=>"Anterodorsal preoptic nucleus", "abbreviation"=>"ADP", "color"=>"FF5547"}, {"id"=>226, "name"=>"Lateral preoptic area", "abbreviation"=>"LPO", "color"=>"F2483B"}, {"id"=>258, "name"=>"Lateral septal nucleus, rostral (rostroventral) part", "abbreviation"=>"LSr", "color"=>"90CBED"}, {"id"=>342, "name"=>"Substantia innominata", "abbreviation"=>"SI", "color"=>"A2B1D8"}, {"id"=>351, "name"=>"Bed nuclei of the stria terminalis", "abbreviation"=>"BST", "color"=>"B3C0DF"}, {"id"=>452, "name"=>"Median preoptic nucleus", "abbreviation"=>"MEPO", "color"=>"FF5547"}, {"id"=>523, "name"=>"Medial preoptic area", "abbreviation"=>"MPO", "color"=>"FF5547"}, {"id"=>564, "name"=>"Medial septal nucleus", "abbreviation"=>"MS", "color"=>"96A7D3"}, {"id"=>596, "name"=>"Diagonal band nucleus", "abbreviation"=>"NDB", "color"=>"96A7D3"}, {"id"=>1109, "name"=>"Parastrial nucleus", "abbreviation"=>"PS", "color"=>"FF5547"}]
## product.id sum num.voxels structure_level_1 structure_level_2
## 7 5 18.51436 478177 Cerebrum Cerebral nuclei
## structure_level_3 tracer_sum tracer_max_distance
## 7 Pallidum 177476.1 8.412521tracer_image <- "/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/stats/data/projection/158257355/projection_density.mnc"
tracer_vol <- mincArray(mincGetVolume(tracer_image))
anatFile <- "/hpf/largeprojects/MICe/edeguzman/misscanada/multi_modal_data/transformed_mnc_files/mt/mtr_average.mnc"
anatVolume <- mincArray(mincGetVolume(anatFile))
mincPlotSliceSeries(anatVolume, tracer_vol, anatLow=0, anatHigh=0.2, low=0.1, high=0.5, symmetric=F, begin=50, end=-40, legend="tracer density", plottitle = "", col=colorRampPalette(c("black", "green"))(100))