multivariate-microarrays

Meta-analysis is a valuable tool for the synthesis of evidence across a wide range study types including high-throughput experiments such as genome-wide association studies (GWAS) and gene expression studies. There are situations though, in which we have multiple outcomes or multiple treatments, in which the multivariate meta-analysis framework which performs a joint modeling of the different quantities of interest may offer important advantages, such as increasing statistical power and allowing performing global tests. Here we apply for the first time these methods to data regarding inflammatory bowel disease (IBD), with its two main forms, Crohn’s disease (CD) and Ulcerative colitis (UC), sharing similar symptoms, but differing in the location and extent of inflammation and in complications.The meta-analysis was performed on 4 datasets containing data for 25,409 unique genes/loci from 251 CD patients, 175 UC patients and 159 healthy individuals. We identified 249 statistically significant differentially expressed genes (DEGs) in CD and 38 in UC at an FDR of 1%. 20 of the DEGs were found be common in both disorders. With the use of the global test derived by the multivariate method we also identified 260 statistically significant DEGs, 53 of which were not found in any of the disorders. Pathway analysis revealed important molecular pathways implicated in the pathogenesis of IBD, such as the JAK/STAT and Interferon-gamma signaling pathways, genes involved in carcinogenesis and others. The multivariate method will be useful in future applications. 




set maxvar 32767

use "test13.dta",clear

set more off

file open meta using mvmeta.txt, write append

file open metan using single_study.txt, write append

file write meta "gene"
file write meta " "
file write meta "b1"
file write meta " "
file write meta "se_b1"
file write meta " "
file write meta "p_b1" 
file write meta " "
file write meta "b2"
file write meta " "
file write meta "se_b2"
file write meta " "
file write meta "p_b2" 
file write meta " "
file write meta "global1"
file write meta " "
file write meta "global2"_n

file write metan "gene"
file write metan " "
file write metan "b1"
file write metan " "
file write metan "se_b1"
file write metan " "
file write metan "z_b1"
file write metan " "
file write metan "p_b1"
file write metan " "
file write metan  "df_b1"
file write metan " "
file write metan "b2"
file write metan " "
file write metan "se_b2"
file write metan " "
file write metan "z_b2"
file write metan " "
file write metan "p_b2"
file write metan " "
file write metan  "df_b2"
file write metan " "
file write metan "w1"
file write metan " "
file write metan "global1"
file write metan " "
file write metan "w2"
file write metan " "
file write metan "global2" _n

  qui sum study
  local x=r(max)

foreach var of varlist a1cf - znf888 {
preserve

     qui gen last=0
     qui gen x0=0
     qui gen x1=0
     qui gen x2=0
     qui gen n0=0
     qui gen n1=0
     qui gen n2=0
     qui gen sd0=0
     qui gen sd1=0
     qui gen sd2=0
     qui gen d1=0
     qui gen d2=0
     qui gen v11=0
     qui gen v22=0
     qui gen v12=0
     qui gen n=0
     qui gen m=0
       qui gen double gmfm=0
     qui gen double gmfm_1=0
     qui gen J=0
     qui gen b1=0
        qui gen b2=0
     qui gen V11=0
     qui gen V22=0
     qui gen V12=0
     qui gen se1=0
     qui gen se2=0


    forvalues i=1/`x' {
   
            qui sum `var' if study==`i' & status==0
            qui replace x0=r(mean) if study==`i'
            qui replace n0=r(N)     if study==`i'
            qui replace sd0=r(sd) if study==`i'   
   
            qui  sum `var' if study==`i' & status==1
            qui replace x1=r(mean) if study==`i'
            qui replace n1=r(N)     if study==`i'
            qui replace sd1=r(sd) if study==`i'
           
            qui sum `var' if study==`i' & status==2
            qui replace x2=r(mean) if study==`i'
            qui replace n2=r(N)     if study==`i'
            qui replace sd2=r(sd) if study==`i'   
       
             qui replace n= (n1+ n2+ n0) if study==`i'
             qui replace m= (n1+ n2+ n0 -3) if study==`i'
             qui replace d1=(x1 - x0)/(sqrt( (    ((n1 - 1)*(sd1^2)) +((n2 - 1)*(sd2^2)) + ((n0 - 1)*(sd0^2)))/(m))) if study==`i'
             qui replace d2=(x2 - x0)/(sqrt( (    ((n1 - 1)*(sd1^2)) +((n2 - 1)*(sd2^2)) + ((n0 - 1)*(sd0^2)))/(m))) if study==`i'
             qui replace v11= (1/n1) + (1/n0) + ((d1^2)/(2*(n))) if study==`i'
             qui replace v22= (1/n2) + (1/n0) + ((d2^2)/(2*(n))) if study==`i'
             qui replace v12= (1/n0) +((d1*d2)/(2*(n))) if study==`i'
             qui replace gmfm = exp((lngamma(m/2))) if study==`i'
             qui replace gmfm_1 = exp((lngamma((m-1)/2))) if study==`i'
             qui replace J = gmfm/((sqrt(m/2))*gmfm_1) if study==`i'
             qui replace b1 = (J*(d1)) if study==`i'             
             qui replace b2 = (J*(d2)) if study==`i'
             qui replace V11 = ((J^2)*(v11)) if study==`i'
             qui replace V22 = ((J^2)*(v22)) if study==`i'
             qui replace V12 = ((J^2)*(v12)) if study==`i'
             qui replace se1=sqrt(V11)
             qui replace se2=sqrt(V22)
       
}

bysort study : replace last=_n==_N
         keep if last==1


   qui summarize x0
   qui return list

if (r(N)==1) {

    keep if x0!=.
   
      qui   metan b1 se1
     local var_es1= r(ES)
     local var_sees1= r(seES)
     local var_z1= r(z)
     local var_p1= r(p_z)
     local var_df1= r(df)
     
       qui metan b2 se2
     local var_es2= r(ES)
     local var_sees2= r(seES)
     local var_z2= r(z)
     local var_p2= r(p_z)
     local var_df2= r(df)
     
    qui mat d=(b1[1],b2[1])
    qui mat list d
    qui mat g=(b1[1]\b2[1])
    qui mat list g
    qui mat v=(V11[1],V12[1]\V12[1],V22[1])
    qui mat list v
    qui mat V=invsym(v)
    qui mat list V
    qui mat W=d*V*g
    qui mat list W
    local var_w1=W[1,1]
    local var_global1=1-chi2(2, W[1,1])
   
    qui mat w2=(b1[1]-b2[1])/ (sqrt( (V11[1]+V22[1])-(2*V12[1])))
    qui mat list w2
    local var_w2=w2[1,1]
    local var_global2=2*(1-normal(abs( w2[1,1] )))
   
   dis `var_es1', `var_sees1', `var_z1', `var_p1', `var_df1', `var_es2', `var_sees2', `var_z2', `var_p2', `var_df2', `var_w1', `var_w2', `var_global1', `var_global2'

        file write  metan "`var'"
        file write  metan " "
        file write  metan "`var_es1'"
        file write  metan " "
        file write  metan "`var_sees1'"
        file write  metan " "
        file write  metan "`var_z1'"
        file write  metan " "
        file write  metan "`var_p1'"
        file write  metan " "
        file write  metan  "`var_df1'"
        file write  metan " "
        file write  metan "`var_es2'"
        file write  metan " "
        file write  metan "`var_sees2'"
        file write  metan " "
        file write  metan "`var_z2'"
        file write  metan " "
        file write  metan "`var_p2'"
        file write  metan " "
        file write  metan  "`var_df2'"
        file write  metan " "
        file write  metan  "`var_w1'"
        file write  metan " "
        file write  metan "`var_global1'"
        file write  metan " "
        file write  metan  "`var_w2'"
        file write  metan " "
        file write  metan "`var_global2'"_n  
    }

    else {
        qui mvmeta b V, vars(b1 b2) mm
          qui matrix v=e(V)
       
           qui matrix a=e(b)

           qui matrix c=e(V)

        ** coef b1
           local var_b1=a[1,1]
       

        ** coef b2
           local var_b2=a[1,2]
       

        ** std.err sto tetragono b1
           local var_st_b1=c[1,1]
          
        ** std.err b1 
        local se_b1=sqrt(c[1,1])
       
        ** std.err sto tetragono b2
           local var_st_b2=c[2,2]
          
         ** std.err b2
          local se_b2=sqrt(c[1,1])
         
         
        **calculation of p-value b1
           local var_p_b1=2*(1-normal(abs(a[1,1]/(sqrt(c[1,1])))))
       
   
        **calculation of p-value b2
           local var_p_b2=2*(1-normal(abs(a[1,2]/(sqrt(c[2,2])))))
       
       
           qui test b1 b2
           qui return list
           local var_g1= r(p)
         
           qui test b1=b2
           qui return list
           local var_g2= r(p)
                   
        dis `var_b1', `se_b1', `var_p_b1', `var_b2', `se_b2', `var_p_b2', `var_g1', `var_g2'

        file write  meta "`var'"
        file write  meta " "
        file write  meta "`var_b1'"
        file write  meta " "
        file write  meta "`se_b1'"
        file write  meta " "
        file write  meta "`var_p_b1'"
        file write  meta " "
        file write  meta "`var_b2'"
        file write  meta " "
        file write  meta  "`se_b2'"
        file write  meta " "
        file write  meta "`var_p_b2'"
        file write  meta " "
        file write  meta "`var_g1'"
        file write  meta " "
        file write  meta "`var_g2'"_n
        }
restore
}

file close metan

file close meta

ċ
multivariate_code_13.do
(7k)
Pantelis Bagos,
Aug 1, 2019, 11:47 AM
ċ
test13.dta
(40k)
Pantelis Bagos,
Aug 1, 2019, 11:46 AM
Comments