Gene network analysis for MyConnectome data

Code available at: https://github.com/poldrack/myconnectome/blob/master/myconnectome/rnaseq/Run_WGCNA.Rmd

This code loads the variance stabilized data (residualized against RIN and the first 3 PCs across subjects) and performs weighted gene coexpression network analysis (WGCNA). See the WGCNA page at UCLA for more info on this technique.

## ==========================================================================
## *
## *  Package WGCNA 1.47 loaded.
## *
## ==========================================================================
## Warning in allowWGCNAThreads: Requested number of threads is higher than number
##  of available processors (or cores). Using too many threads may degrade code performance. It is recommended that the number of threads is no more than number
##  of available processors.
## 
## Allowing multi-threading with up to 2 threads.
## pickSoftThreshold: will use block size 7370.
##  pickSoftThreshold: calculating connectivity for given powers...
##    ..working on genes 1 through 7370 of 13844
##    ..working on genes 7371 through 13844 of 13844
##    Power SFT.R.sq  slope truncated.R.sq  mean.k. median.k. max.k.
## 1      1   0.0359  0.675          0.932 2670.000  2.65e+03 4110.0
## 2      2   0.2730 -1.250          0.912  806.000  7.65e+02 1860.0
## 3      3   0.5840 -1.770          0.949  309.000  2.71e+02 1010.0
## 4      4   0.6970 -2.030          0.956  138.000  1.09e+02  617.0
## 5      5   0.7520 -2.090          0.963   69.000  4.80e+01  403.0
## 6      6   0.7870 -2.070          0.963   37.600  2.26e+01  277.0
## 7      7   0.8380 -1.950          0.970   21.900  1.12e+01  197.0
## 8      8   0.9000 -1.850          0.985   13.500  5.82e+00  148.0
## 9      9   0.9210 -1.860          0.993    8.690  3.14e+00  122.0
## 10    10   0.9370 -1.840          0.996    5.820  1.75e+00  102.0
## 11    12   0.9570 -1.770          0.993    2.880  5.85e-01   76.0
## 12    14   0.9690 -1.680          0.989    1.580  2.12e-01   58.5
## 13    16   0.9630 -1.620          0.979    0.938  8.22e-02   46.3
## 14    18   0.9790 -1.560          0.994    0.593  3.37e-02   37.9
## 15    20   0.9830 -1.500          0.996    0.394  1.46e-02   31.5

Plot of network topology results to examine the scale-free topology metric:

plot of chunk plotPickThresh

Based on these results, we choose power=8 to obtain fit of r2~0.9

##  Calculating module eigengenes block-wise from all genes
##    Flagging genes and samples with too many missing values...
##     ..step 1
##  ....pre-clustering genes to determine blocks..
##    Projective K-means:
##    ..k-means clustering..
##    ..merging smaller clusters...
## Block sizes:
## gBlocks
##    1    2    3    4 
## 3994 3991 3723 2136 
##  ..Working on block 1 .
##     TOM calculation: adjacency..
##     ..will use 2 parallel threads.
##      Fraction of slow calculations: 0.000000
##     ..connectivity..
##     ..matrix multiplication..
##     ..normalization..
##     ..done.
##  ....clustering..
##  ....detecting modules..
##  ....calculating module eigengenes..
##  ....checking modules for statistical meaningfulness..
##      ..removing 1 genes from module 1 because their KME is too low.
##      ..removing 1 genes from module 13 because their KME is too low.
##  ..Working on block 2 .
##     TOM calculation: adjacency..
##     ..will use 2 parallel threads.
##      Fraction of slow calculations: 0.000000
##     ..connectivity..
##     ..matrix multiplication..
##     ..normalization..
##     ..done.
##  ....clustering..
##  ....detecting modules..
##  ....calculating module eigengenes..
##  ....checking modules for statistical meaningfulness..
##      ..removing 6 genes from module 1 because their KME is too low.
##  ..Working on block 3 .
##     TOM calculation: adjacency..
##     ..will use 2 parallel threads.
##      Fraction of slow calculations: 0.000000
##     ..connectivity..
##     ..matrix multiplication..
##     ..normalization..
##     ..done.
##  ....clustering..
##  ....detecting modules..
##  ....calculating module eigengenes..
##  ....checking modules for statistical meaningfulness..
##      ..removing 1 genes from module 4 because their KME is too low.
##      ..removing 1 genes from module 10 because their KME is too low.
##  ..Working on block 4 .
##     TOM calculation: adjacency..
##     ..will use 2 parallel threads.
##      Fraction of slow calculations: 0.000000
##     ..connectivity..
##     ..matrix multiplication..
##     ..normalization..
##     ..done.
##  ....clustering..
##  ....detecting modules..
##  ....calculating module eigengenes..
##  ....checking modules for statistical meaningfulness..
##   ..reassigning 2 genes from module 1 to modules with higher KME.
##   ..reassigning 1 genes from module 2 to modules with higher KME.
##   ..reassigning 2 genes from module 3 to modules with higher KME.
##   ..reassigning 1 genes from module 4 to modules with higher KME.
##   ..reassigning 1 genes from module 5 to modules with higher KME.
##   ..reassigning 1 genes from module 6 to modules with higher KME.
##   ..reassigning 1 genes from module 8 to modules with higher KME.
##   ..reassigning 1 genes from module 9 to modules with higher KME.
##   ..reassigning 1 genes from module 10 to modules with higher KME.
##   ..reassigning 1 genes from module 12 to modules with higher KME.
##   ..reassigning 1 genes from module 15 to modules with higher KME.
##   ..reassigning 1 genes from module 20 to modules with higher KME.
##   ..reassigning 33 genes from module 40 to modules with higher KME.
##   ..reassigning 15 genes from module 41 to modules with higher KME.
##   ..reassigning 7 genes from module 42 to modules with higher KME.
##   ..reassigning 13 genes from module 43 to modules with higher KME.
##   ..reassigning 2 genes from module 44 to modules with higher KME.
##   ..reassigning 2 genes from module 46 to modules with higher KME.
##   ..reassigning 1 genes from module 48 to modules with higher KME.
##   ..reassigning 1 genes from module 49 to modules with higher KME.
##   ..reassigning 2 genes from module 50 to modules with higher KME.
##   ..reassigning 1 genes from module 51 to modules with higher KME.
##   ..reassigning 1 genes from module 52 to modules with higher KME.
##   ..reassigning 1 genes from module 54 to modules with higher KME.
##   ..reassigning 2 genes from module 70 to modules with higher KME.
##   ..reassigning 2 genes from module 72 to modules with higher KME.
##   ..reassigning 1 genes from module 79 to modules with higher KME.
##  ..merging modules that are too close..
##      mergeCloseModules: Merging modules whose distance is less than 0.15
##        Calculating new MEs...

Plot WGCNA results

plot of chunk plotResults

## png 
##   2

get hub genes using module membership (KME)

V1 V2 V3 V4 V5 V6
kME1 ANKRD28 PRKG1 NCKAP1 CD109 C7orf41
kME2 TSPO SPSB3 OXLD1 CYBA HMHA1
kME3 RBBP6 RAD50 SAP130 SMARCC2 ZBTB4
kME4 BPTF PRRC2C PHF3 NCOR1 SPEN
kME5 BRD8 KIAA0040 AHCYL1 UTP14C HMG20A
kME6 PLXNC1 CLEC7A HCK NCF2 TNFAIP2
kME7 CREBBP CHD2 SPTAN1 KAT6A EP400
kME8 C14orf1 PTBP1 ZNF490 TRIM4 SSR1
kME9 HNRNPM CERS5 GGA2 TADA2B LAS1L
kME10 ZNF609 NASP EIF3A IK SAFB
kME11 EXOC1 ST8SIA4 MARCH1 KCNJ2 ZFYVE16
kME12 THAP5 PTPRCAP ERO1L CHURC1 ZNF624
kME13 SNHG8 CCDC109B RPS19BP1 EIF3E ZNF277
kME14 CCDC25 SH3KBP1 RTF1 FLOT2 SUPT16H
kME15 FAM196B ANGPT2 TTN TENM1 A2M
kME16 RGS14 PPP1R12C PCED1A LIMD2 MTA1
kME17 JMJD1C TCEB3 CCDC88A CD68 TSC1
kME18 FKBP15 ZCCHC6 OSCAR SIRPB1 TRIM21
kME19 ZHX2 PACS1 TTC3 CYFIP2 DOCK9
kME20 HNMT AIF1 CSTA TNFSF13B CPVL
kME21 FLII HCLS1 PI4KB STK10 EHMT1
kME22 DERL1 SEC23B DPY19L4 UNC50 MAGT1
kME23 METTL9 BZW1 UQCRQ TMEM167B PSMA1
kME24 CNTRL KTN1 EEA1 TPR CAST
kME25 SEPHS2 HIF1AN HNRNPA2B1 RAB27A APBB1IP
kME26 GAA CORO1B SH3TC1 SLC16A3 HK3
kME27 BMP2K FGD4 PLXDC2 CPEB4 ZNF117
kME28 LRRN3 FAIM3 DYRK2 ATP8B2 GPR183
kME29 PFDN1 HADHA URI1 BLMH ZCCHC7
kME30 GLB1 ZDHHC7 PSEN1 ARSD TOR1B
kME31 SH2D1B AKR1C3 KLRD1 PRR5L MLC1
kME32 ATP5A1 WDR61 IER3IP1 GHITM ATP5B
kME33 BTN2A1 PRPF38B ZDHHC5 ELF4 IL17RA
kME34 RBM3 CHMP1B VOPP1 ATXN10 SPG21
kME35 UBR4 ANKRD36B TNRC6B ITPR2 HECTD4
kME36 FAM107B RALA SAP30L TMED10 UBE2D2
kME37 PARP9 RSAD2 IFIT3 MX1 IFIT2
kME38 CEP78 NCAM1 S1PR5 ATP2B4 COLQ
kME39 SCAF4 WASF2 HNRNPD TARDBP MDC1
kME40 H6PD MTF1 BLOC1S3 ZDHHC5 IPO7
kME41 ABLIM1 IDS SYNRG CDC25B CBX5
kME42 CMIP MITF ASAP1 MAML3 CYB5R3
kME43 UCP2 CLSTN1 PDCL3 ARL4C HNRNPA0
kME44 VTI1A STRN4 DHX8 RAPGEF2 REL
kME45 C12orf49 PHF5A CRLS1 MRPL49 AGPAT6
kME46 VMA21 HMGN4 TMEM8A NUDT21 TMEM248
kME47 LRRC25 LPCAT2 MS4A7 SCARB2 TMTC1
kME48 TBX21 GNLY KIR3DL1 USP28 CTSW
kME49 RRAGA SCAMP2 COQ5 SS18L2 SERBP1
kME50 PPP2R3C PSMA6 PPP2R4 PSMD8 MRPS22
kME51 UNC50 TMEM30A RNFT1 MIR17HG STARD3NL
kME52 BCL2A1 RNF149 PNPLA8 RELT ATG16L2
kME53 LRRK2 DMXL2 WDFY3 AQP9 GLUL
kME54 P4HB ACAA2 SUMF1 UQCRC1 PDIA3
kME55 MS4A3 HDC CCR3 CLC CPA3
kME56 ABRACL TNFRSF14 PPP1R12C C6orf211 HNRNPF
kME57 KIAA0430 XPO7 GSK3B IDS HIPK1
kME58 ZBTB11 SENP7 ZNF92 ESCO1 RAB11FIP2
kME59 CEP63 MIA3 DNAJC7 RAB5C TBCA
kME60 COMMD8 MRPS35 FAU STXBP3 MRPL13
kME61 MDH1 RAN IARS2 OSTC DDOST
kME62 FBXL5 PGK1 OAZ2 NDUFS1 LY86
kME63 LTF OLFM4 CEACAM8 CRISP3 BPI