setwd("/dartfs-hpc/rc/home/m/f0052zm/szhao_lab_share/jihyun/diff")

foldernames <- list.dirs(recursive=FALSE)

foldernames <- foldernames[ grepl("20220316", foldernames) ]

foldernames

1. './BLCA_20220316'
2. './BRCA_20220316'
3. './CESC_20220316'
4. './CHOL_20220316'
5. './ESCA_20220316'
6. './GBM_20220316'
7. './HNSC_20220316'
8. './KICH_20220316'
9. './KIRC_20220316'
10. './KIRP_20220316'
11. './LIHC_20220316'
12. './LUAD_20220316'
13. './LUSC_20220316'
14. './PAAD_20220316'
15. './PRAD_20220316'
16. './SKCM_20220316'
17. './TGCT_20220316'
18. './UCEC_20220316'
19. './UCS_20220316'

Significant genes identified by diffdriver

19 tumor types were analyzed, each run for 46 traits. Use Fdr 0.1 as cut off

outdf <- NULL
for (i in 1: length(foldernames)){
  tumorName <- gsub("\\./", "", gsub("\\_20220316","",foldernames[i]))
  fileNames <- dir(file.path(getwd(),foldernames[i]),pattern = "resdd.txt", full.names = F) 
    for (j in 1:length(fileNames)){
    data <- read.table(file.path(foldernames[i],fileNames[j]), comment="", header = TRUE)
    data <- data[data$dd.fdr < 0.1, ]
    if (dim(data)[1]!= 0){
        data$trait <- gsub("\\_res.*","", fileNames[j])
        data$gene <- rownames(data)
        data$tumor <- tumorName
        outdf <- rbind(data, outdf)
        }
    }
}
rownames(outdf) <- NULL
outdf <- outdf[with(outdf, order(tumor, dd.fdr)), ]
outdf <- outdf[, c(9,8,7,1:6)]

Add parameter estimates

outdf2 <- cbind(outdf, data.frame("alpha.null" = NA, "alpha0" = NA, "alpha1" = NA))
for (i in 1:nrow(outdf2)){
    #for (i in 1){
        load(paste0(outdf2[i, "tumor"], "_20220316/",  outdf2[i, "trait"], "_resdd.", outdf2[i, "gene"], ".Rd"))
        par <- diffdriver::ddmodel(mutmtx, canno$Phenotype, bmrmtx, fe[,1])
        outdf2[i, c("alpha.null", "alpha0", "alpha1")] <- c(par$null.param, par$alt.param)
     }

outdf2

	tumor	gene	trait	dd	mlr	mlr.v2	dd.fdr	mlr.fdr	mlr.v2.fdr	alpha.null	alpha0	alpha1
32	BLCA	FMN2	SmokingCessation	7.425869e-05	0.0006166599	0.0004852809	0.004678297	0.03884957	0.03057270	0.262500	27.802814	-835.80883
34	BLCA	ARID1B	DrinksPerWeek	4.235875e-04	0.0208948295	0.0217697748	0.026686013	0.48318215	0.49090826	0.887500	737.223764	1812.75205
31	BLCA	EP300	SmokingCessation	1.385923e-03	0.0391951293	0.0299120078	0.043656587	0.61732329	0.47111412	2.700000	274.044797	-534.03808
28	BLCA	PHF3	ukbEUR_VARICOSE_VEINS	7.362218e-04	0.0814185435	0.1326985317	0.046381976	0.59184271	0.70483140	-0.462500	-991.712252	344.23066
33	BLCA	SF1	EA_NG_2018_excluding_23andMe	7.492196e-04	0.0172975203	0.0201645437	0.047200832	0.44911674	0.54473980	0.446875	-1268.026471	-835.51528
35	BLCA	FMN2	DrinksPerWeek	2.009992e-03	0.0100022931	0.0105232748	0.063314755	0.48318215	0.49090826	0.262500	736.761944	1783.85048
29	BLCA	MKI67	SWB_NG_2015	1.290483e-03	0.0274749483	0.0355558091	0.081300459	0.47528141	0.59501867	-1.075000	-184.536061	1342.90636
36	BLCA	HSP90AA1	BIP_Cell_2018	1.350536e-03	0.1122305068	0.1794930847	0.085083747	0.71725528	0.93424350	1.000000	-334.322712	276.58283
30	BLCA	ARID1B	SmokingCessation	4.289904e-03	0.2857414967	0.3021672147	0.090087992	0.94988805	0.97494509	0.887500	13.657937	522.11792
26	HNSC	THSD7A	ukbEUR_PSYCHIATRIC	3.713252e-04	0.0031069903	0.0008982485	0.012253730	0.10253068	0.02964220	-1.037500	-75.606495	1206.52550
24	HNSC	FAT1	ukbEUR_PSYCHIATRIC	2.626886e-03	0.3856815163	0.4371952893	0.043343622	0.79220879	0.80882955	25.500000	175.746976	-541.54579
23	HNSC	EPHA2	ukbEUR_VARICOSE_VEINS	2.389289e-03	0.0031548965	0.0014798485	0.078846533	0.10411159	0.04883500	1.181250	-833.734965	297.73784
27	HNSC	SMARCA4	ukbEUR_ALLERGIC_RHINITIS	2.783025e-03	0.1491054680	0.1298596654	0.091839819	0.61580552	0.61219557	0.562500	321.519477	-447.50115
25	HNSC	RASA1	ukbEUR_PSYCHIATRIC	9.003776e-03	0.7425501030	0.7516081161	0.099041536	0.87514834	0.91863214	1.196875	173.313675	824.14750
22	KIRC	ATM	IMMU	8.112777e-03	0.2257837733	0.2483776339	0.081127765	0.53587523	0.58213486	4.150000	41.999920	412.08851
18	LIHC	ARID1A	ukbEUR_OBESITY	3.901772e-03	0.1034694733	0.0931384191	0.050723030	0.60095231	0.58453082	0.156250	501.266724	297.69201
19	LIHC	ARID1A	ukbEUR_BMI	5.660597e-03	0.1369189521	0.1285657704	0.073587758	0.60418585	0.55711834	0.156250	921.636060	573.96375
20	LIHC	ARID1A	SWB_NG_2015	9.031156e-03	0.0267209791	0.0283458154	0.074498598	0.31857859	0.36849560	0.156250	-116.007523	1320.90844
21	LIHC	COL11A1	SWB_NG_2015	1.146132e-02	0.0490120912	0.0593100587	0.074498598	0.31857859	0.38551538	-0.956250	-163.687454	-1579.80320
16	LUAD	MGA	BMD_2018	1.714333e-03	0.0015073537	0.0009960654	0.028819881	0.02863972	0.01892524	-0.078125	386.461637	436.09637
17	LUAD	RB1	BMD_2018	3.033672e-03	0.1794517589	0.1789524305	0.028819881	0.54514936	0.58855214	0.243750	358.729732	473.50387
14	LUAD	DZIP1L	ukbEUR_CASES	2.170081e-03	0.0108658180	0.0127687754	0.041231546	0.11342601	0.12130337	0.865625	4.182155	705.56767
13	LUAD	RB1	ukbEUR_TI_cojo	2.624951e-03	0.2917339348	0.2791322239	0.049874067	0.46191206	0.61820304	0.243750	-360.937163	-981.84793
15	LUAD	SETD2	SWB_NG_2015	4.497941e-03	0.0560942966	0.0667918080	0.085460878	0.48588359	0.48814607	3.700000	91.829827	1068.50944
12	LUSC	ADAMTS12	ukbEUR_CASES	5.834716e-03	0.0118340775	0.0083106424	0.058347164	0.11834078	0.08310642	38.300000	-18.716907	-517.44985
11	PRAD	ATM	ukbEUR_HERNIA_ABDOMINOPELVIC	7.072951e-03	0.1314269383	0.1352235289	0.099021315	0.52776395	0.51851314	1.125000	366.455047	-564.06793
8	SKCM	HNF4G	IBD_NG2017.hap3	4.454584e-04	0.0195460393	0.0130187890	0.008909168	0.21458073	0.20533229	1.675000	110.131365	-178.36080
10	SKCM	HNF4G	BIP_Cell_2018	5.092831e-04	0.0102335756	0.0103167574	0.010185662	0.12603606	0.13558585	1.675000	392.272450	-222.77599
5	SKCM	SORCS3	ukbEUR_CI_cojo	8.312359e-04	0.3803239603	0.1389332665	0.016624717	0.69149811	0.48607939	-1.062500	-28.211284	-88.18981
2	SKCM	HNF4G	ukbEUR_TI_cojo	1.603427e-03	0.6132987439	0.5337082699	0.032068540	0.97823085	0.95085321	1.675000	445.982925	983.11689
9	SKCM	TACR3	Breast_Cancer	2.315632e-03	0.1511012484	0.2144934872	0.046312640	0.93354579	0.88919101	1.575000	130.017749	250.59202
3	SKCM	GUCA1C	ukbEUR_Height	3.172305e-03	0.0192193267	0.0124527862	0.063446105	0.38438653	0.24905572	-0.262500	799.914457	-426.00289
4	SKCM	GUCA1C	ukbEUR_DYSLIPID	3.376700e-03	0.0380163428	0.0520147899	0.067533993	0.38016343	0.26690634	-0.262500	128.909225	-283.87926
7	SKCM	CDH9	RA_NG2010.hap3	4.412549e-03	0.0177085564	0.0061885981	0.088250980	0.17708556	0.06188598	3.050000	375.187072	-238.15489
6	SKCM	TACR3	UC_NG_2015	4.987284e-03	0.0101304365	0.0149340136	0.099745682	0.08508095	0.09956009	1.575000	-1432.698813	-203.78165
1	UCEC	EP300	ukbEUR_DYSLIPID	2.089489e-03	0.0161016066	0.0284157264	0.083579551	0.32203213	0.55214288	0.218750	115.919507	-384.57865

Add plot for each gene

options(repr.plot.width =7, repr.plot.height =3)
for (i in 1:nrow(outdf)){
    load(paste0(outdf[i, "tumor"], "_20220316/",  outdf[i, "trait"], "_resdd.", outdf[i, "gene"], ".Rd"))
    diffdriver::plot_mut(mutmtx, canno, bmrmtx, ganno, main = paste(outdf[i,c("tumor","gene", "trait")], collapse = "-"))
}

options(repr.plot.width =10, repr.plot.height =6)
labels <- apply(outdf2[, c("gene", "trait")], 1, paste0, collapse = "-")
par(mar=c(10,4,4,4))
a <- barplot(-log10(t(as.matrix(outdf2[, c("dd", "mlr.v2")]))), col= c("salmon","darkgreen"), beside = T, ylab = "-log10(p)", xaxt ='n')
text(a[1,], par("usr")[3], labels = apply(outdf2[, c("gene", "trait")], 1, paste0, collapse = "-"), srt = 45, adj = c(1.1,1.1), xpd = TRUE, cex=0.6)