1b_DESeq2_SM

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("DESeq2")

library(Matrix)
library(data.table)
library(dplyr)
library(DESeq2)

RNGkind("L'Ecuyer-CMRG")

data.dir = "/scratch/as5916/Projects/Rutgers_asthma_smc_scSeq/Data"

# ------------------------------------------------------------------------
# read in cell information
# ------------------------------------------------------------------------

cell_info = fread(file.path(data.dir, "meta_lung.tsv"), 
                  stringsAsFactors=TRUE)
dim(cell_info)
cell_info[1:2,]

read in processed count data as a single matrix

# ------------------------------------------------------------------------
# read in count data of one celltype
# ------------------------------------------------------------------------

dat = readRDS(file.path(data.dir, sprintf("ct_mtx/combinedCountMatrix.rds"))) #, grp
dim(dat)

[1]  25146 173551

class(dat)

[1] "dgCMatrix"
attr(,"package")
[1] "Matrix"

dat[1:5,1:4]

5 x 4 sparse Matrix of class "dgCMatrix"
                AAACCAAAGAACCTAT-1_1 AAACCAAAGATTGCAT-1_1 AAACCAAAGATTGCGC-1_1 AAACCAAAGGCTACTA-1_1
ENSG00000238009                    .                    .                    .                    .
ENSG00000241860                    .                    .                    .                    .
ENSG00000290385                    .                    .                    .                    .
ENSG00000291215                    .                    .                    .                    .
LINC01409                          .                    .                    .                    .

gc()

            used   (Mb) gc trigger   (Mb)   max used    (Mb)
Ncells   7520060  401.7   11612667  620.2   11612667   620.2
Vcells 879478856 6709.9 1269052289 9682.2 3738572678 28523.1

————————————————————————

subset cell information

————————————————————————

table(colnames(dat) %in% cell_info$cell)

  TRUE 
173551 

meta = cell_info[match(colnames(dat), cell_info$cell),]
dim(meta) #173551     19

[1] 173551     19

meta[1:2,]

#names(meta_cell)[11:12] = c("PMI", "RIN")
#names(meta_cell)[15:16] = c("mitoPercent", "riboPercent")
#dim(meta_cell)
#meta_cell[1:2,]

#summary(meta)
#meta_cell$age = scale(meta_cell$age)
#meta_cell$PMI = scale(meta_cell$PMI)
#meta_cell$RIN = scale(meta_cell$RIN)
#meta_cell$Capbatch = droplevels(meta_cell$Capbatch)
#meta_cell$Seqbatch = droplevels(meta_cell$Seqbatch)
#meta_cell$individual = as.factor(meta_cell$individual)

#meta_cell$age = scale(meta$Aage)

meta$subject_ID = as.factor(meta$subject_ID)

adjust certain column names in meta2kp to match the requirement

of ideas and for the ease of later processing

meta2kp <- meta

gc()

            used   (Mb) gc trigger   (Mb)   max used    (Mb)
Ncells   7519744  401.6   11612667  620.2   11612667   620.2
Vcells 879353253 6709.0 1269052289 9682.2 3738572678 28523.1

original 1b_DESeq2 does not have this section - why??

# adjust certain column names in meta2kp to match the requirement 
# of ideas and for the ease of later processing

# colnames_meta2kp = names(meta2kp)
colnames_meta2kp = names(meta2kp)
names(meta2kp)[which(colnames_meta2kp=="cell")] = "cell_id"

#names(meta2kp)[which(colnames_meta2kp=="donor")] = "individual"
names(meta2kp)[which(colnames_meta2kp=="subject_ID")] = "individual"

names(meta2kp)[which(colnames_meta2kp=="Disease")] = "diagnosis"

names(meta2kp)[which(colnames_meta2kp=="Sex")] = "sex"

names(meta2kp)[which(colnames_meta2kp=="Age")] = "age"

#'age',
meta_ind = distinct(meta2kp[,c('individual', 'diagnosis',  'age', 'sex')])
table(meta_ind$diagnosis, meta_ind$sex)

        
         Female Male
  Asthma      1    3
  Normal      1    3

————————————————————————

collect count data

————————————————————————

trec1 = matrix(NA, nrow=nrow(dat), ncol=nrow(meta_ind))
colnames(trec1) = meta_ind$subject_ID

rownames(trec1) = rownames(dat)
dim(trec1)

[1] 25146     8

trec1[1:2,1:3]

                [,1] [,2] [,3]
ENSG00000238009   NA   NA   NA
ENSG00000241860   NA   NA   NA

trec1 = matrix(NA, nrow=nrow(dat), ncol=nrow(meta_ind)) colnames(trec1) = meta_ind$subject_ID

rownames(trec1) = rownames(dat) dim(trec1) [1] 25146 8 trec1[1:2,1:3] [,1] [,2] [,3] ENSG00000238009 NA NA NA ENSG00000241860 NA NA NA

#dat$subject_ID
for(i in 1:ncol(trec1)){
  wi = which(meta2kp$individual == meta_ind$individual[i])
  #trec1[,i] = rowSums(dat1[,wi])
  trec1[,i] = rowSums(dat[,wi])
}

dim(trec1)

[1] 25146     8

trec1[1:2,1:3]

                [,1] [,2] [,3]
ENSG00000238009   51   58   49
ENSG00000241860  194  237  176

for(i in 1:ncol(trec1)){ + wi = which(meta2kp\(individual == meta_ind\)individual[i]) + #trec1[,i] = rowSums(dat1[,wi]) + trec1[,i] = rowSums(dat[,wi]) + }

dim(trec1) [1] 25146 8 trec1[1:2,1:3] [,1] [,2] [,3] ENSG00000238009 51 58 49 ENSG00000241860 194 237 176

————————————————————————

run DESeq2

————————————————————————

colData = meta_ind
colnames(colData)[2] = 'diagnosis'
# keep the age column being numeric
for(i in 1:(ncol(colData)-1)){
  if(is.character(colData[[i]])){
    colData[[i]] = as.factor(colData[[i]])
  }
}
dim(colData)

[1] 8 4

colData[1:2,]

NA

colData = meta_ind colnames(colData)[2] = ‘diagnosis’ # keep the age column being numeric for(i in 1:(ncol(colData)-1)){ + if(is.character(colData[[i]])){ + colData[[i]] = as.factor(colData[[i]]) + } + } dim(colData) [1] 8 3 colData[1:2,] individual diagnosis sex 1: N04091813 Normal Female 2: A510101913 Asthma Female

summary(colData)

      individual  diagnosis      age           sex   
 A502061813:1    Asthma:4   Min.   :11.0   Female:2  
 A504262313:1    Normal:4   1st Qu.:14.0   Male  :6  
 A507182213:1               Median :17.0             
 A510101913:1               Mean   :17.5             
 N011118413:1               3rd Qu.:22.0             
 N04091813 :1               Max.   :23.0             
 (Other)   :2                                        

summary(colData) individual diagnosis age sex
A502061813:1 Asthma:4 Min. :11.0 Female:2
A504262313:1 Normal:4 1st Qu.:14.0 Male :6
A507182213:1 Median :17.0
A510101913:1 Mean :17.5
N011118413:1 3rd Qu.:22.0
N04091813 :1 Max. :23.0
(Other) :2

#colData$diagnosis = factor(colData$diagnosis, levels=c("mild", "severe"))
colData$diagnosis = factor(colData$diagnosis, levels=c("Normal", "Asthma"))

————————————————————————

add to colData donor level covariates to adjust for

————————————————————————

total read depth across all cells under each individual

total_rd = c()

for (i in c(1:dim(meta_ind)[1])){
  wi = which(meta2kp$individual == meta_ind$individual[i])
  #wi = which(meta2kp$donor == meta_ind$donor[i])
  #total_rd = c(total_rd, sum(apply(dat1[,wi], 2, sum)))
  total_rd = c(total_rd, sum(apply(dat[,wi], 2, sum)))
}

colData$totalrd = total_rd

for (i in c(1:dim(meta_ind)[1])){ + wi = which(meta2kp\(individual == meta_ind\)individual[i]) + #wi = which(meta2kp\(donor == meta_ind\)donor[i]) + #total_rd = c(total_rd, sum(apply(dat1[,wi], 2, sum))) + total_rd = c(total_rd, sum(apply(dat[,wi], 2, sum))) + } Warning in asMethod(object) : sparse->dense coercion: allocating vector of size 4.7 GiB Warning in asMethod(object) : sparse->dense coercion: allocating vector of size 4.3 GiB Warning in asMethod(object) : sparse->dense coercion: allocating vector of size 4.1 GiB Warning in asMethod(object) : sparse->dense coercion: allocating vector of size 4.2 GiB Warning in asMethod(object) : sparse->dense coercion: allocating vector of size 3.0 GiB Warning in asMethod(object) : sparse->dense coercion: allocating vector of size 3.3 GiB Warning in asMethod(object) : sparse->dense coercion: allocating vector of size 4.3 GiB Warning in asMethod(object) : sparse->dense coercion: allocating vector of size 4.6 GiB

colData$totalrd = total_rd

first, null model

rm(dat)

dd0 = DESeqDataSetFromMatrix(countData = trec1, 
                             colData = colData,
                             design = ~ diagnosis)
dd0  = DESeq(dd0)

dd0 = DESeqDataSetFromMatrix(countData = trec1, + colData = colData, + design = ~ diagnosis) Warning in S4Vectors:::anyMissing(runValue(x_seqnames)) : ‘S4Vectors:::anyMissing()’ is deprecated. Use ‘anyNA()’ instead. See help(“Deprecated”) Warning in S4Vectors:::anyMissing(runValue(strand(x))) : ‘S4Vectors:::anyMissing()’ is deprecated. Use ‘anyNA()’ instead. See help(“Deprecated”) converting counts to integer mode dd0 = DESeq(dd0) estimating size factors estimating dispersions gene-wise dispersion estimates mean-dispersion relationship final dispersion estimates fitting model and testing

dd0 = LargeDESeqDataset

res0 = results(dd0)
dim(res0)

[1] 25146     6

head(res0)

log2 fold change (MLE): diagnosis Asthma vs Normal 
Wald test p-value: diagnosis Asthma vs Normal 
DataFrame with 6 rows and 6 columns
                 baseMean log2FoldChange     lfcSE       stat    pvalue      padj
                <numeric>      <numeric> <numeric>  <numeric> <numeric> <numeric>
ENSG00000238009   44.5001      0.1716105  0.349439  0.4911034  0.623353   0.99998
ENSG00000241860  181.5093      0.0318759  0.281603  0.1131947  0.909876   0.99998
ENSG00000290385  161.3508     -0.0211281  0.319833 -0.0660597  0.947330   0.99998
ENSG00000291215  493.8922      0.1484186  0.250274  0.5930242  0.553165   0.99998
LINC01409       2259.0088      0.0204456  0.140824  0.1451854  0.884565   0.99998
ENSG00000290784  302.8115     -0.0585569  0.131504 -0.4452862  0.656113   0.99998

summary(res0)

out of 25146 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 14, 0.056%
LFC < 0 (down)     : 19, 0.076%
outliers [1]       : 130, 0.52%
low counts [2]     : 0, 0%
(mean count < 1)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

res0 = results(dd0) dim(res0) [1] 25146 6 head(res0) log2 fold change (MLE): diagnosis Asthma vs Normal Wald test p-value: diagnosis Asthma vs Normal DataFrame with 6 rows and 6 columns baseMean log2FoldChange lfcSE stat pvalue padj ENSG00000238009 44.5001 0.1716105 0.349439 0.4911034 0.623353 0.99998 ENSG00000241860 181.5093 0.0318759 0.281603 0.1131947 0.909876 0.99998 ENSG00000290385 161.3508 -0.0211281 0.319833 -0.0660597 0.947330 0.99998 ENSG00000291215 493.8922 0.1484186 0.250274 0.5930242 0.553165 0.99998 LINC01409 2259.0088 0.0204456 0.140824 0.1451854 0.884565 0.99998 ENSG00000290784 302.8115 -0.0585569 0.131504 -0.4452862 0.656113 0.99998 summary(res0)

out of 25146 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 14, 0.056% LFC < 0 (down) : 19, 0.076% outliers [1] : 130, 0.52% low counts [2] : 0, 0% (mean count < 1) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results

save

saveRDS(dd0, file = "/scratch/as5916/Projects/Rutgers_asthma_smc_scSeq/res/deseq2/modeldd0.rds")
saveRDS(res0, file = "/scratch/as5916/Projects/Rutgers_asthma_smc_scSeq/res/deseq2/modeldd0_res0.rds")

nm0 = resultsNames(dd0)
nm0

nm0 = resultsNames(dd0) nm0 [1] “Intercept” “diagnosis_Asthma_vs_Normal”

nm0 = nm0[-1]

nm0 [1] “diagnosis_Asthma_vs_Normal”

pvals0 = matrix(NA, nrow=nrow(trec1), ncol=length(nm0))

for(k in 1:length(nm0)){
  rk = results(dd0, name=nm0[k])
  pvals0[,k] = rk$pvalue
}

colnames(pvals0) = nm0
dim(pvals0)
head(pvals0)
summary(pvals0)

colnames(pvals0) = nm0 dim(pvals0) [1] 25146 1 head(pvals0) diagnosis_Asthma_vs_Normal [1,] 0.6233533 [2,] 0.9098762 [3,] 0.9473303 [4,] 0.5531649 [5,] 0.8845645 [6,] 0.6561129 summary(pvals0) diagnosis_Asthma_vs_Normal Min. :0.0000
1st Qu.:0.3147
Median :0.5727
Mean :0.5520
3rd Qu.:0.7959
Max. :1.0000
NA’s :130

second, only include donor level total read depth as covariate

what is log(totalrd) = total read depth

dd1 = DESeqDataSetFromMatrix(countData = trec1, 
                             colData = colData,
                             design = ~ log(totalrd) + diagnosis)
dd1 = DESeq(dd1)

dd1 = DESeqDataSetFromMatrix(countData = trec1, + colData = colData, + design = ~ log(totalrd) + diagnosis) Warning in S4Vectors:::anyMissing(runValue(x_seqnames)) : ‘S4Vectors:::anyMissing()’ is deprecated. Use ‘anyNA()’ instead. See help(“Deprecated”) Warning in S4Vectors:::anyMissing(runValue(strand(x))) : ‘S4Vectors:::anyMissing()’ is deprecated. Use ‘anyNA()’ instead. See help(“Deprecated”) converting counts to integer mode the design formula contains one or more numeric variables with integer values, specifying a model with increasing fold change for higher values. did you mean for this to be a factor? if so, first convert this variable to a factor using the factor() function the design formula contains one or more numeric variables that have mean or standard deviation larger than 5 (an arbitrary threshold to trigger this message). Including numeric variables with large mean can induce collinearity with the intercept. Users should center and scale numeric variables in the design to improve GLM convergence. dd1 = DESeq(dd1) estimating size factors estimating dispersions gene-wise dispersion estimates mean-dispersion relationship final dispersion estimates fitting model and testing 45 rows did not converge in beta, labelled in mcols(object)$betaConv. Use larger maxit argument with nbinomWaldTest

res1 = results(dd1)
dim(res1)

[1] 25146     6

head(res1)

log2 fold change (MLE): diagnosis Asthma vs Normal 
Wald test p-value: diagnosis Asthma vs Normal 
DataFrame with 6 rows and 6 columns
                 baseMean log2FoldChange     lfcSE       stat    pvalue      padj
                <numeric>      <numeric> <numeric>  <numeric> <numeric> <numeric>
ENSG00000238009   44.5001      0.1676752 0.3947612  0.4247510  0.671018  0.999995
ENSG00000241860  181.5093      0.0552185 0.2745308  0.2011378  0.840591  0.999995
ENSG00000290385  161.3508      0.0107464 0.3420452  0.0314181  0.974936  0.999995
ENSG00000291215  493.8922      0.0799207 0.2263560  0.3530751  0.724032  0.999995
LINC01409       2259.0088     -0.0122597 0.0776207 -0.1579435  0.874501  0.999995
ENSG00000290784  302.8115     -0.0628565 0.1485842 -0.4230361  0.672269  0.999995

summary(res1)

out of 25146 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 677, 2.7%
LFC < 0 (down)     : 617, 2.5%
outliers [1]       : 0, 0%
low counts [2]     : 0, 0%
(mean count < 1)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

res1 = results(dd1) dim(res1) [1] 25146 6 head(res1) log2 fold change (MLE): diagnosis Asthma vs Normal Wald test p-value: diagnosis Asthma vs Normal DataFrame with 6 rows and 6 columns baseMean log2FoldChange lfcSE stat pvalue padj ENSG00000238009 44.5001 0.1676752 0.3947612 0.4247510 0.671018 0.999995 ENSG00000241860 181.5093 0.0552185 0.2745308 0.2011378 0.840591 0.999995 ENSG00000290385 161.3508 0.0107464 0.3420452 0.0314181 0.974936 0.999995 ENSG00000291215 493.8922 0.0799207 0.2263560 0.3530751 0.724032 0.999995 LINC01409 2259.0088 -0.0122597 0.0776207 -0.1579435 0.874501 0.999995 ENSG00000290784 302.8115 -0.0628565 0.1485842 -0.4230361 0.672269 0.999995 summary(res1)

out of 25146 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 677, 2.7% LFC < 0 (down) : 617, 2.5% outliers [1] : 0, 0% low counts [2] : 0, 0% (mean count < 1) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results

save

saveRDS(dd1, file = "/scratch/as5916/Projects/Rutgers_asthma_smc_scSeq/res/deseq2/modeldd1.rds")
saveRDS(res1, file = "/scratch/as5916/Projects/Rutgers_asthma_smc_scSeq/res/deseq2/modeldd1_res1.rds")

nm1 = resultsNames(dd1)
nm1

[1] "Intercept"                  "log.totalrd."               "diagnosis_Asthma_vs_Normal"

nm1 = nm1[-1]
nm1

[1] "log.totalrd."               "diagnosis_Asthma_vs_Normal"

“log.totalrd.” “diagnosis_Asthma_vs_Normal”

pvals1 = matrix(NA, nrow=nrow(trec1), ncol=length(nm1))

for(k in 1:length(nm1)){
  rk = results(dd1, name=nm1[k])
  pvals1[,k] = rk$pvalue
}

colnames(pvals1) = nm1
dim(pvals1)

[1] 25146     2

head(pvals1)

     log.totalrd. diagnosis_Asthma_vs_Normal
[1,] 0.8532099534                  0.6710182
[2,] 0.2215338620                  0.8405908
[3,] 0.5474593215                  0.9749361
[4,] 0.1151434580                  0.7240321
[5,] 0.0002114418                  0.8745013
[6,] 0.7701444422                  0.6722689

summary(pvals1)

  log.totalrd.       diagnosis_Asthma_vs_Normal
 Min.   :3.000e-08   Min.   :0.0000            
 1st Qu.:2.740e-01   1st Qu.:0.3046            
 Median :5.582e-01   Median :0.5690            
 Mean   :5.351e-01   Mean   :0.5395            
 3rd Qu.:8.169e-01   3rd Qu.:0.7951            
 Max.   :9.999e-01   Max.   :1.0000            

colnames(pvals1) = nm1 dim(pvals1) [1] 25146 2 head(pvals1) log.totalrd. diagnosis_Asthma_vs_Normal [1,] 0.8532099534 0.6710182 [2,] 0.2215338620 0.8405908 [3,] 0.5474593215 0.9749361 [4,] 0.1151434580 0.7240321 [5,] 0.0002114418 0.8745013 [6,] 0.7701444422 0.6722689 summary(pvals1) log.totalrd. diagnosis_Asthma_vs_Normal Min. :3.000e-08 Min. :0.0000
1st Qu.:2.740e-01 1st Qu.:0.3046
Median :5.582e-01 Median :0.5690
Mean :5.351e-01 Mean :0.5395
3rd Qu.:8.169e-01 3rd Qu.:0.7951
Max. :9.999e-01 Max. :1.0000

third, include sex and age as covariates

Sex + Age + Dx

dd2 = DESeqDataSetFromMatrix(countData = trec1, 
                             colData = colData,
                             design = ~ sex + age + diagnosis)
dd2 = DESeq(dd2)

dd2 = DESeqDataSetFromMatrix(countData = trec1, + colData = colData, + design = ~ sex + age + diagnosis) Warning in S4Vectors:::anyMissing(runValue(x_seqnames)) : ‘S4Vectors:::anyMissing()’ is deprecated. Use ‘anyNA()’ instead. See help(“Deprecated”) Warning in S4Vectors:::anyMissing(runValue(strand(x))) : ‘S4Vectors:::anyMissing()’ is deprecated. Use ‘anyNA()’ instead. See help(“Deprecated”) converting counts to integer mode the design formula contains one or more numeric variables with integer values, specifying a model with increasing fold change for higher values. did you mean for this to be a factor? if so, first convert this variable to a factor using the factor() function the design formula contains one or more numeric variables that have mean or standard deviation larger than 5 (an arbitrary threshold to trigger this message). Including numeric variables with large mean can induce collinearity with the intercept. Users should center and scale numeric variables in the design to improve GLM convergence. dd2 = DESeq(dd2) estimating size factors estimating dispersions gene-wise dispersion estimates mean-dispersion relationship final dispersion estimates fitting model and testing

res2 = results(dd2)
dim(res2)

[1] 25146     6

head(res2)

log2 fold change (MLE): diagnosis Asthma vs Normal 
Wald test p-value: diagnosis Asthma vs Normal 
DataFrame with 6 rows and 6 columns
                 baseMean log2FoldChange     lfcSE       stat    pvalue      padj
                <numeric>      <numeric> <numeric>  <numeric> <numeric> <numeric>
ENSG00000238009   44.5001      0.1594950  0.424258  0.3759392  0.706962  0.999977
ENSG00000241860  181.5093      0.0123845  0.320962  0.0385856  0.969221  0.999977
ENSG00000290385  161.3508     -0.0213933  0.297474 -0.0719167  0.942668  0.999977
ENSG00000291215  493.8922      0.0985977  0.238604  0.4132272  0.679440  0.999977
LINC01409       2259.0088      0.0151563  0.121460  0.1247840  0.900695  0.999977
ENSG00000290784  302.8115     -0.0594934  0.146352 -0.4065089  0.684369  0.999977

summary(res2)

out of 25146 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 47, 0.19%
LFC < 0 (down)     : 53, 0.21%
outliers [1]       : 0, 0%
low counts [2]     : 1463, 5.8%
(mean count < 4)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

res2 = results(dd2) dim(res2) [1] 25146 6 head(res2) log2 fold change (MLE): diagnosis Asthma vs Normal Wald test p-value: diagnosis Asthma vs Normal DataFrame with 6 rows and 6 columns baseMean log2FoldChange lfcSE stat pvalue padj ENSG00000238009 44.5001 0.1594950 0.424258 0.3759392 0.706962 0.999977 ENSG00000241860 181.5093 0.0123845 0.320962 0.0385856 0.969221 0.999977 ENSG00000290385 161.3508 -0.0213933 0.297474 -0.0719167 0.942668 0.999977 ENSG00000291215 493.8922 0.0985977 0.238604 0.4132272 0.679440 0.999977 LINC01409 2259.0088 0.0151563 0.121460 0.1247840 0.900695 0.999977 ENSG00000290784 302.8115 -0.0594934 0.146352 -0.4065089 0.684369 0.999977 summary(res2)

out of 25146 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 47, 0.19% LFC < 0 (down) : 53, 0.21% outliers [1] : 0, 0% low counts [2] : 1463, 5.8% (mean count < 4) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results

save

saveRDS(dd2, file = "/scratch/as5916/Projects/Rutgers_asthma_smc_scSeq/res/deseq2/modeldd2.rds")
saveRDS(res2, file = "/scratch/as5916/Projects/Rutgers_asthma_smc_scSeq/res/deseq2/modeldd2_res2.rds")

nm2 = resultsNames(dd2)
nm2

[1] "Intercept"                  "sex_Male_vs_Female"         "age"                        "diagnosis_Asthma_vs_Normal"

nm2 = resultsNames(dd2) nm2 [1] “Intercept” “sex_Male_vs_Female” “age” “diagnosis_Asthma_vs_Normal”

nm2 = nm2[-1]
nm2

[1] "sex_Male_vs_Female"         "age"                        "diagnosis_Asthma_vs_Normal"

nm2 [1] “sex_Male_vs_Female” “age” “diagnosis_Asthma_vs_Normal”

pvals2 = matrix(NA, nrow=nrow(trec1), ncol=length(nm2))

for(k in 1:length(nm2)){
  rk = results(dd2, name=nm2[k])
  pvals2[,k] = rk$pvalue
}

dim(pvals2)

[1] 25146     3

head(pvals2)

           [,1]      [,2]      [,3]
[1,] 0.54767163 0.7094586 0.7069621
[2,] 0.53425963 0.6699251 0.9692208
[3,] 0.14058974 0.5334659 0.9426682
[4,] 0.23248508 0.1610069 0.6794402
[5,] 0.05961753 0.2183867 0.9006945
[6,] 0.53383192 0.2004347 0.6843688

summary(pvals2)

       V1               V2               V3        
 Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
 1st Qu.:0.1387   1st Qu.:0.2782   1st Qu.:0.3209  
 Median :0.3938   Median :0.5191   Median :0.5694  
 Mean   :0.4228   Mean   :0.5175   Mean   :0.5457  
 3rd Qu.:0.6856   3rd Qu.:0.7708   3rd Qu.:0.7865  
 Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  

dim(pvals2) [1] 25146 3 head(pvals2) [,1] [,2] [,3] [1,] 0.54767163 0.7094586 0.7069621 [2,] 0.53425963 0.6699251 0.9692208 [3,] 0.14058974 0.5334659 0.9426682 [4,] 0.23248508 0.1610069 0.6794402 [5,] 0.05961753 0.2183867 0.9006945 [6,] 0.53383192 0.2004347 0.6843688 summary(pvals2) V1 V2 V3
Min. :0.0000 Min. :0.0000 Min. :0.0000
1st Qu.:0.1387 1st Qu.:0.2782 1st Qu.:0.3209
Median :0.3938 Median :0.5191 Median :0.5694
Mean :0.4228 Mean :0.5175 Mean :0.5457
3rd Qu.:0.6856 3rd Qu.:0.7708 3rd Qu.:0.7865
Max. :1.0000 Max. :1.0000 Max. :1.0000

fourth, include both donor level total read depth, sex and age

Sex + Age + Dx

dd3 = DESeqDataSetFromMatrix(countData = trec1, 
                             colData = colData,
                             design = ~ log(totalrd) + sex + age + diagnosis)
dd3 = DESeq(dd3)

dd3 = DESeqDataSetFromMatrix(countData = trec1, + colData = colData, + design = ~ log(totalrd) + sex + age + diagnosis) Warning in S4Vectors:::anyMissing(runValue(x_seqnames)) : ‘S4Vectors:::anyMissing()’ is deprecated. Use ‘anyNA()’ instead. See help(“Deprecated”) Warning in S4Vectors:::anyMissing(runValue(strand(x))) : ‘S4Vectors:::anyMissing()’ is deprecated. Use ‘anyNA()’ instead. See help(“Deprecated”) converting counts to integer mode the design formula contains one or more numeric variables with integer values, specifying a model with increasing fold change for higher values. did you mean for this to be a factor? if so, first convert this variable to a factor using the factor() function the design formula contains one or more numeric variables that have mean or standard deviation larger than 5 (an arbitrary threshold to trigger this message). Including numeric variables with large mean can induce collinearity with the intercept. Users should center and scale numeric variables in the design to improve GLM convergence. dd3 = DESeq(dd3) estimating size factors estimating dispersions gene-wise dispersion estimates mean-dispersion relationship – note: fitType=‘parametric’, but the dispersion trend was not well captured by the function: y = a/x + b, and a local regression fit was automatically substituted. specify fitType=‘local’ or ‘mean’ to avoid this message next time. final dispersion estimates fitting model and testing 4 rows did not converge in beta, labelled in mcols(object)$betaConv. Use larger maxit argument with nbinomWaldTest

res3 = results(dd3)
dim(res3)

[1] 25146     6

head(res3)

log2 fold change (MLE): diagnosis Asthma vs Normal 
Wald test p-value: diagnosis Asthma vs Normal 
DataFrame with 6 rows and 6 columns
                 baseMean log2FoldChange     lfcSE       stat    pvalue      padj
                <numeric>      <numeric> <numeric>  <numeric> <numeric> <numeric>
ENSG00000238009   44.5001     0.15774965 0.3873647  0.4072381  0.683833  0.999999
ENSG00000241860  181.5093     0.03811529 0.2555879  0.1491279  0.881453  0.999999
ENSG00000290385  161.3508     0.00350413 0.3019385  0.0116054  0.990740  0.999999
ENSG00000291215  493.8922     0.08072926 0.2470384  0.3267883  0.743828  0.999999
LINC01409       2259.0088    -0.00636910 0.0478966 -0.1329760  0.894212  0.999999
ENSG00000290784  302.8115    -0.04893732 0.1037474 -0.4716970  0.637143  0.999999

summary(res3)

out of 25146 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 419, 1.7%
LFC < 0 (down)     : 398, 1.6%
outliers [1]       : 0, 0%
low counts [2]     : 0, 0%
(mean count < 1)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

res3 = results(dd3) dim(res3) [1] 25146 6 head(res3) log2 fold change (MLE): diagnosis Asthma vs Normal Wald test p-value: diagnosis Asthma vs Normal DataFrame with 6 rows and 6 columns baseMean log2FoldChange lfcSE stat pvalue padj ENSG00000238009 44.5001 0.15774965 0.3873647 0.4072381 0.683833 0.999999 ENSG00000241860 181.5093 0.03811529 0.2555879 0.1491279 0.881453 0.999999 ENSG00000290385 161.3508 0.00350413 0.3019385 0.0116054 0.990740 0.999999 ENSG00000291215 493.8922 0.08072926 0.2470384 0.3267883 0.743828 0.999999 LINC01409 2259.0088 -0.00636910 0.0478966 -0.1329760 0.894212 0.999999 ENSG00000290784 302.8115 -0.04893732 0.1037474 -0.4716970 0.637143 0.999999 summary(res3)

out of 25146 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 419, 1.7% LFC < 0 (down) : 398, 1.6% outliers [1] : 0, 0% low counts [2] : 0, 0% (mean count < 1) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results

save

saveRDS(dd3, file = "/scratch/as5916/Projects/Rutgers_asthma_smc_scSeq/res/deseq2/modeldd3.rds")
saveRDS(res3, file = "/scratch/as5916/Projects/Rutgers_asthma_smc_scSeq/res/deseq2/modeldd3_res3.rds")

nm3 = resultsNames(dd3)
nm3

[1] "Intercept"                  "log.totalrd."               "sex_Male_vs_Female"         "age"                        "diagnosis_Asthma_vs_Normal"

nm3 [1] “Intercept” “log.totalrd.” “sex_Male_vs_Female” “age” “diagnosis_Asthma_vs_Normal”

nm3 = nm3[-1]
nm3

[1] "log.totalrd."               "sex_Male_vs_Female"         "age"                        "diagnosis_Asthma_vs_Normal"

nm3 [1] “log.totalrd.” “sex_Male_vs_Female” “age” “diagnosis_Asthma_vs_Normal”

pvals3 = matrix(NA, nrow=nrow(trec1), ncol=length(nm3))

for(k in 1:length(nm3)){
  rk = results(dd3, name=nm3[k])
  pvals3[,k] = rk$pvalue
}

colnames(pvals3) = nm3
dim(pvals3)

[1] 25146     4

head(pvals3)

     log.totalrd. sex_Male_vs_Female        age diagnosis_Asthma_vs_Normal
[1,] 0.6596686991         0.42048228 0.92466368                  0.6838331
[2,] 0.0991147453         0.84091133 0.15748019                  0.8814527
[3,] 0.7224578509         0.29331661 0.48120531                  0.9907404
[4,] 0.5556609110         0.49513670 0.41864739                  0.7438280
[5,] 0.0001679168         0.04330512 0.51021470                  0.8942123
[6,] 0.4236755584         0.23992254 0.05113285                  0.6371431

summary(pvals3)

  log.totalrd.    sex_Male_vs_Female      age         diagnosis_Asthma_vs_Normal
 Min.   :0.0000   Min.   :0.0000     Min.   :0.0000   Min.   :0.0000            
 1st Qu.:0.3875   1st Qu.:0.2110     1st Qu.:0.2836   1st Qu.:0.3173            
 Median :0.6666   Median :0.5251     Median :0.5868   Median :0.6198            
 Mean   :0.6087   Mean   :0.5108     Mean   :0.5502   Mean   :0.5682            
 3rd Qu.:0.8753   3rd Qu.:0.8151     3rd Qu.:0.8282   3rd Qu.:0.8368            
 Max.   :0.9999   Max.   :1.0000     Max.   :1.0000   Max.   :1.0000            

colnames(pvals3) = nm3 dim(pvals3) [1] 25146 4 head(pvals3) log.totalrd. sex_Male_vs_Female age diagnosis_Asthma_vs_Normal [1,] 0.6596686991 0.42048228 0.92466368 0.6838331 [2,] 0.0991147453 0.84091133 0.15748019 0.8814527 [3,] 0.7224578509 0.29331661 0.48120531 0.9907404 [4,] 0.5556609110 0.49513670 0.41864739 0.7438280 [5,] 0.0001679168 0.04330512 0.51021470 0.8942123 [6,] 0.4236755584 0.23992254 0.05113285 0.6371431 summary(pvals3) log.totalrd. sex_Male_vs_Female age diagnosis_Asthma_vs_Normal Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
1st Qu.:0.3875 1st Qu.:0.2110 1st Qu.:0.2836 1st Qu.:0.3173
Median :0.6666 Median :0.5251 Median :0.5868 Median :0.6198
Mean :0.6087 Mean :0.5108 Mean :0.5502 Mean :0.5682
3rd Qu.:0.8753 3rd Qu.:0.8151 3rd Qu.:0.8282 3rd Qu.:0.8368
Max. :0.9999 Max. :1.0000 Max. :1.0000 Max. :1.0000

results

pdf(sprintf("/scratch/as5916/Projects/Rutgers_asthma_smc_scSeq/res/deseq2/1b_DESeq2_pval_hists.pdf"), 
    width=18, height=3.5)
par(mfrow=c(1,4), bty="n", mar=c(5,4,2,1))
k = length(nm0)
hist(pvals0[,k], main="null", xlab="p-value", breaks=50)
k = length(nm1)
hist(pvals1[,k], main="log(totalrd)", xlab="p-value", breaks=50)

k = length(nm2)
hist(pvals2[,k], main="sex+age", xlab="p-value", breaks=50)
k = length(nm3)
hist(pvals3[,k], main="log(totalrd)+sex+age", xlab="p-value", breaks=50)

dev.off()

null device 
          1 

par(mfrow=c(1,4), bty="n", mar=c(5,4,2,1))
k = length(nm0)
hist(pvals0[,k], main="null", xlab="p-value", breaks=50)
k = length(nm1)
hist(pvals1[,k], main="log(totalrd)", xlab="p-value", breaks=50)

k = length(nm2)
hist(pvals2[,k], main="sex+age", xlab="p-value", breaks=50)
k = length(nm3)
hist(pvals3[,k], main="log(totalrd)+sex+age", xlab="p-value", breaks=50)

pdf(sprintf("/scratch/as5916/Projects/Rutgers_asthma_smc_scSeq/res/deseq2/1b_DESeq2_pval_hists_full.pdf"), 
    width=9, height=7)
par(mfrow=c(2,2), bty="n", mar=c(5,4,2,1))
for(k in 1:length(nm3)){
  hist(pvals3[,k], main=nm3[k], xlab="p-value", breaks=50)
}

dev.off()

null device 
          1 

par(mfrow=c(2,2), bty="n", mar=c(5,4,2,1))
for(k in 1:length(nm3)){
  hist(pvals3[,k], main=nm3[k], xlab="p-value", breaks=50)
}

NA
NA

————————————————————————

summarize p-value distribution

————————————————————————

pdf(sprintf("/scratch/as5916/Projects/Rutgers_asthma_smc_scSeq/res/deseq2/1b_DESeq2_compare_pval.pdf"), 
    width=9, height=3)
par(mfrow=c(1,3), bty="n", mar=c(5,4,1,1))

plot(-log10(res0$pvalue), -log10(res1$pvalue), main="-log10(p-value)", 
     xlab="null", ylab="log(totalrd)", 
     pch=20, cex=0.2, col=rgb(0.8, 0.1, 0.1, 0.5))
abline(0, 1, col="darkblue")

plot(-log10(res0$pvalue), -log10(res2$pvalue), main="-log10(p-value)", 
     xlab="null", ylab="sex+age", 
     pch=20, cex=0.2, col=rgb(0.8, 0.1, 0.1, 0.5))

abline(0, 1, col="darkblue")

plot(-log10(res0$pvalue), -log10(res3$pvalue), main="-log10(p-value)", 
     xlab="null", ylab="log(totalrd)+sex+age", 
     pch=20, cex=0.2, col=rgb(0.8, 0.1, 0.1, 0.5))
abline(0, 1, col="darkblue")

dev.off()

null device 
          1 

par(mfrow=c(1,3), bty="n", mar=c(5,4,1,1))

plot(-log10(res0$pvalue), -log10(res1$pvalue), main="-log10(p-value)", 
     xlab="null", ylab="log(totalrd)", 
     pch=20, cex=0.2, col=rgb(0.8, 0.1, 0.1, 0.5))
abline(0, 1, col="darkblue")

plot(-log10(res0$pvalue), -log10(res2$pvalue), main="-log10(p-value)", 
     xlab="null", ylab="sex+age", 
     pch=20, cex=0.2, col=rgb(0.8, 0.1, 0.1, 0.5))
abline(0, 1, col="darkblue")

plot(-log10(res0$pvalue), -log10(res3$pvalue), main="-log10(p-value)", 
     xlab="null", ylab="log(totalrd)+sex+age", 
     pch=20, cex=0.2, col=rgb(0.8, 0.1, 0.1, 0.5))
abline(0, 1, col="darkblue")