import pandas as pd
import numpy as np
import modas.multiprocess as mp
from rpy2.robjects.packages import importr
from rpy2.robjects import pandas2ri
from rpy2.rinterface_lib.embedded import RRuntimeError
import rpy2.robjects as robjects
from collections import Counter
import modas.gwas_cmd as gc
import pyranges as pr
from yattag import Doc, indent
import subprocess
import shutil
import warnings
import glob
import os
import re
pandas2ri.activate()
data_table = importr('data.table')
base = importr('base')
robjects.r('options(datatable.showProgress = FALSE)')
warnings.filterwarnings("ignore")
utils_path = subprocess.check_output('locate modas/utils', shell=True, text=True, encoding='utf-8')
#utils_path = '/'.join(re.search('\n(.*site-packages.*)\n', utils_path).group(1).split('/')[:-1])
utils_path = re.search('\n(.*site-packages.*)\n', utils_path).group(1)
if not utils_path.endswith('utils'):
utils_path = '/'.join(utils_path.split('/')[:-1])
# def gwas(phe, geno, num_threads, phe_fn):
# geno_prefix = geno.split('/')[-1]
# related_matrix_cmd = 'gemma.linux -bfile {0}.link -gk 1 -o {1}'.format(geno_prefix,geno_prefix)
# gwas_cmd = 'gemma.linux -bfile {0}.link -k output/{0}.cXX.txt -lmm -n {1} -o {2}'
# fam = pd.read_csv(geno+'.fam', sep=r'\s+', header=None)
# fam[5] = 1
# fam = pd.merge(fam, phe, left_on=0, right_index=True, how='left')
# fam.to_csv(geno_prefix+'.link.fam', sep='\t', na_rep='NA', header=None, index=False)
# if os.path.exists(geno_prefix+'.link.bed'):
# os.remove(geno_prefix+'.link.bed')
# if os.path.exists(geno_prefix+'.link.bim'):
# os.remove(geno_prefix+'.link.bim')
# os.symlink(geno+'.bed', geno_prefix+'.link.bed')
# os.symlink(geno+'.bim', geno_prefix+'.link.bim')
# values = list()
# for _, p in enumerate(phe.columns):
# p = p.replace('/', '.')
# values.append((gwas_cmd.format(*[geno_prefix, _ + 2, '.'.join(phe_fn.split('/')[-1].split('.')[:-1])+ '_' + str(p)]),))
# s = mp.run(related_matrix_cmd)
# if s != 0:
# return None
# else:
# s = mp.parallel(mp.run, values, num_threads)
# os.remove(geno_prefix+'.link.bed')
# os.remove(geno_prefix+'.link.bim')
# os.remove(geno_prefix+'.link.fam')
# return s
def gwas(phe, geno, num_threads, phe_fn, gwas_model):
software, model = gwas_model.split('_')
geno_prefix = geno.split('/')[-1]
phe_fn = '.'.join(phe_fn.split('/')[-1].split('.')[:-1])
if software == 'gemma' and model == 'MLM':
geno_prefix = geno.split('/')[-1]
related_matrix_cmd = utils_path+'/gemma.linux -bfile {0}.link -gk 1 -o {1}'.format(geno_prefix,geno_prefix)
fam = pd.read_csv(geno + '.fam', sep=r'\s+', header=None)
fam[5] = 1
fam = pd.merge(fam, phe, left_on=0, right_index=True, how='left')
fam.to_csv(geno_prefix+'.link.fam', sep='\t', na_rep='NA', header=None, index=False)
if software != 'GAPIT' and gwas_model != 'gemma_LM':
if os.path.exists(geno_prefix+'.link.bed'):
os.remove(geno_prefix+'.link.bed')
if os.path.exists(geno_prefix+'.link.bim'):
os.remove(geno_prefix+'.link.bim')
os.symlink(geno+'.bed', geno_prefix+'.link.bed')
os.symlink(geno+'.bim', geno_prefix+'.link.bim')
if software == 'rMVP':
if os.path.exists(geno_prefix+'.link.fam'):
os.remove(geno_prefix+'.link.fam')
os.symlink(geno + '.fam', geno_prefix + '.link.fam')
fam = pd.read_csv(geno + '.fam', sep=r'\s+', header=None)
phe = phe.reindex(fam[0])
if software == 'gemma' and model == "LM":
if os.path.exists('./gemma_lm_geno'):
shutil.rmtree('./gemma_lm_geno')
os.mkdir('./gemma_lm_geno')
for p in phe.columns:
p = p.replace('m/z', 'm.z')
os.symlink('../' + geno + '.bed', './gemma_lm_geno/' + geno_prefix + '_' + p + '.link.bed')
os.symlink('../' + geno + '.bim', './gemma_lm_geno/' + geno_prefix + '_' + p + '.link.bim')
fam = pd.read_csv(geno + '.fam', sep=r'\s+', header=None)
fam = pd.merge(fam.iloc[:, :5], phe.loc[:, p.replace('m.z', 'm/z')].to_frame(), left_on=0, right_index=True, how='left')
fam.to_csv('./gemma_lm_geno/' + geno_prefix + '_' + p + '.link.fam', sep='\t', na_rep='NA', header=None, index=False)
values = list()
for _, p in enumerate(phe.columns):
p = p.replace('m/z', 'm.z')
if gwas_model == 'gemma_LM':
values.append(((gc.gemma_cmd(model, './gemma_lm_geno/' + geno_prefix + '_' + p + '.link', None, None, '_'.join([phe_fn, gwas_model, p]))),))
if gwas_model == 'gemma_MLM':
values.append(((gc.gemma_cmd(model, geno_prefix + '.link', geno_prefix, _ + 2, '_'.join([phe_fn, gwas_model, p]))),))
if software == 'rMVP':
if not os.path.exists('./output'):
os.mkdir('./output')
omics_phe = phe.loc[:, p.replace('m.z', 'm/z')].to_frame().reset_index()
omics_phe.columns = ['Taxa', '_'.join([phe_fn, gwas_model, p])]
values.append((model, geno_prefix+'.link', geno_prefix+'.link', omics_phe, 1, './output'))
if gwas_model == 'gemma_MLM':
s = mp.run(related_matrix_cmd)
if s != 0:
return None
if software == 'gemma':
s = mp.parallel(mp.run, values, num_threads)
if software == 'rMVP':
s1 = mp.parallel(gc.rmvp, (values[0],), 1)
s = mp.parallel(gc.rmvp, values[1:], num_threads)
s = s1 + s
if software == 'GAPIT':
if not os.path.exists('./output'):
os.mkdir('./output')
phe.columns = ['_'.join([phe_fn, gwas_model, p]) for p in phe.columns]
phe = phe.reset_index()
phe.columns = ['Taxa'] + list(phe.columns[1:])
geno = os.path.abspath(geno)
os.chdir('./output')
s = gc.gapit(model, geno, phe, utils_path)
os.chdir('../')
if gwas_model == 'gemma_MLM' or software == 'rMVP':
os.remove(geno_prefix+'.link.bed')
os.remove(geno_prefix+'.link.bim')
os.remove(geno_prefix+'.link.fam')
if gwas_model == 'gemma_LM':
shutil.rmtree('./gemma_lm_geno')
return s
def gwas_plot(res, p, prefix, t, software):
try:
base.sink('/dev/null')
w = data_table.fread(res, data_table=base.getOption("datatable.fread.datatable", False))
if software == 'gemma':
w_subset = w.loc[w.p_wald <= float(p), :]
m = w_subset[['rs', 'chr', 'ps', 'p_wald']]
q = w[['rs', 'chr', 'ps', 'p_wald']]
if software == 'rMVP':
w_subset = w.loc[w[w.columns[-1]] <= float(p), :]
m = w_subset.iloc[:, [0, 1, 2, -1]]
q = w.iloc[:, [0, 1, 2, -1]]
if software == 'GAPIT':
w_subset = w.loc[w['P.value'] <= float(p), :]
m = w_subset[['SNP', 'Chromosome', 'Position', 'P.value']]
q = w[['SNP', 'Chromosome', 'Position', 'P.value']]
m.columns = ['SNP', 'Chromosome', 'Position', prefix]
q.columns = ['SNP', 'Chromosome', 'Position', prefix]
thresholdi = robjects.FloatVector([1.0 / w.shape[0], 1e-6, 1e-5])
lim = -np.log10(min(m[prefix])) + 2
#w_subset = base.subset(w, np.array(w.rx2('p_wald')) <= float(p))
#m = w_subset.rx(robjects.StrVector(['rs', 'chr', 'ps', 'p_wald']))
#q = w.rx(robjects.StrVector(['rs', 'chr', 'ps', 'p_wald']))
# m.names = ['SNP', 'Chromosome', 'Position', prefix]
# q.names = ['SNP', 'Chromosome', 'Position', prefix]
#thresholdi = robjects.FloatVector([1.0/w.nrow, 1e-6, 1e-5])
#lim = -np.log10(min(np.array(w_subset.rx2('p_wald'))))+2
#base.sink('/dev/null')
robjects.r('source("'+utils_path+'/CMplot.r")')
CMplot = robjects.r['CMplot']
CMplot(m, plot_type='m', col=robjects.StrVector(["grey30", "grey60"]), ylim=robjects.FloatVector([2, lim]), threshold=thresholdi,
cex=robjects.FloatVector([0.5, 0.5, 0.5]), signal_cex=robjects.FloatVector([0.5, 0.5, 0.5]),
threshold_col=robjects.StrVector(['red', 'green', 'blue']), chr_den_col=robjects.rinterface.NULL, amplify=True,
signal_pch = robjects.IntVector([19, 19, 19]), dpi=300,
signal_col=robjects.StrVector(['red', 'green', 'blue']), multracks=False, LOG10=True, file=t)
CMplot(q, plot_type='q', col='grey30', threshold=thresholdi[0],
signal_cex=robjects.FloatVector([0.5, 0.5, 0.5]), signal_pch=robjects.IntVector([19, 19, 19]),
conf_int_col='gray', signal_col='red', multracks=False, LOG10=True, file=t, dpi=300)
base.sink()
except RRuntimeError:
return 0
except ValueError:
return 0
else:
return 1
def gwas_plot_parallel(phe, p, threads, t, phe_fn, gwas_model):
software, model = gwas_model.split('_')
values = list()
for i in phe.columns:
i = str(i).replace('m/z', 'm.z')
if software == 'gemma':
gwas_fn = 'output/' + '_'.join(['.'.join(phe_fn.split('/')[-1].split('.')[:-1]), gwas_model, str(i)]) + '.assoc.txt'
if software == 'rMVP':
gwas_fn = 'output/' + '_'.join(['.'.join(phe_fn.split('/')[-1].split('.')[:-1]), gwas_model, str(i)]) + '.'.join(['.'+model, 'csv'])
if software == 'GAPIT':
gwas_fn = 'output/' + '.'.join(['GAPIT', model, '_'.join(['.'.join(phe_fn.split('/')[-1].split('.')[:-1]), gwas_model, str(i)]), 'GWAS', 'Results', 'csv'])
values.append((gwas_fn, p, '.'.join(phe_fn.split('/')[-1].split('.')[:-1]) + '_' + str(i), t, software))
s = mp.parallel(gwas_plot, values, threads)
return s
def boxplot(phe, g, qtl):
robjects.r('''box_plot <- function(d, phe, rs, level){
library(ggplot2)
library(ggsignif)
d <- d[d$haplotype!=1,]
d[d$haplotype==0,'haplotype'] <- level[1]
d[d$haplotype==2,'haplotype'] <- level[2]
d$haplotype <- factor(d$haplotype,levels = level)
b <- as.numeric(formatC(max(d[,1],na.rm=T)*1.2/4,format = 'e',digits = 1))
p <- ggplot(data = d,aes_string(x='haplotype',y=names(d)[1],fill='haplotype'))+
theme_bw()+
theme(legend.title = element_blank(),
legend.background = element_blank(),
legend.key = element_rect(colour = NA, fill = NA),
legend.text = element_text(size = 4),
legend.key.size = unit(3,'mm'),
legend.position = 'none',
plot.title = element_text(hjust=0.5,size = 6),
plot.margin=unit(c(0.3,0.3,0,0),'cm'),
panel.grid = element_blank(),
axis.line = element_line(colour = 'black',size=0.4),
axis.text = element_text(size = 6,color = 'black'),
axis.ticks.length=unit(.1, 'cm'))+
stat_boxplot(geom = 'errorbar', width = 0.2,size=0.1)+
geom_boxplot(lwd=0.2,width=0.5,outlier.size = 0.2)+
geom_signif(comparisons = list(level),map_signif_level = F,
test= t.test, size=0.2 ,textsize=2, y_position = max(d[,1],na.rm=T)*1.1)+
#xlab('')+ylab('')+ggtitle(paste(phe,rs,sep='_'))+
xlab('')+ylab('')+ggtitle('')+
scale_y_continuous(breaks=seq(0,4*b,by=b),labels = function(x) formatC(x, format = 'e',digits = 1), limits = c(0, max(d[,1], na.rm=T)*1.2))+
scale_fill_manual(values=c('#E3FFE2', 'forest green'))
ggsave(paste(phe,'_',rs,'_','boxplot','.jpg',sep=''),plot=p,device='jpg',width=3.5,height=4.3,units = 'cm')
}''')
robjects.r['options'](warn=-1)
base.sink('/dev/null')
box_plot = robjects.r['box_plot']
g = g.where(g.snp.isin(qtl.SNP), drop=True)
allele = pd.DataFrame([g.a0, g.a1], index=['a0', 'a1'], columns=g.snp.values)
g = pd.DataFrame(g.values, index=g.sample, columns=g.snp.values)
ril = g.index.intersection(phe.index)
g = g.reindex(ril)
phe = phe.reindex(ril)
for index, row in qtl.iterrows():
if row['phe_name'] not in phe.columns:
continue
d = pd.concat([phe[row['phe_name']], g[row['SNP']]], axis=1)
d.columns = ['trait.' + d.columns[0].replace('-', '.'), 'haplotype']
level = robjects.StrVector([allele[row['SNP']]['a1'].values*2, allele[row['SNP']]['a0'].values*2])
box_plot(d, row['phe_name'], row['SNP'], level)
base.sink()
def multi_trait_plot(phe, gwas_dir, qtl, phe_fn, prefix, t, gwas_model):
software, model = gwas_model.split('_')
bk = pd.DataFrame()
for i in phe.columns:
i = i.replace('/', '.')
# fn = gwas_dir + '/' + '.'.join(phe_fn.split('/')[-1].split('.')[:-1]) + '_' + str(i)+'.assoc.txt'
if software == 'gemma':
gwas_fn = gwas_dir + '_'.join(['.'.join(phe_fn.split('/')[-1].split('.')[:-1]), gwas_model, str(i)]) + '.assoc.txt'
d = pd.read_csv(gwas_fn, sep='\t')
d = d[['rs', 'chr', 'ps', 'p_wald']]
if software == 'rMVP':
gwas_fn = gwas_dir + '_'.join(['.'.join(phe_fn.split('/')[-1].split('.')[:-1]), gwas_model, str(i)]) + '.'.join(['.' + model, 'csv'])
d = pd.read_csv(gwas_fn)
d = d.iloc[:, [0, 1, 2, -1]]
d.columns = ['rs', 'chr', 'ps', 'p_wald']
if software == 'GAPIT':
gwas_fn = gwas_dir + '.'.join(['GAPIT', model, '_'.join(['.'.join(phe_fn.split('/')[-1].split('.')[:-1]), gwas_model, str(i)]), 'GWAS', 'Results', 'csv'])
d = pd.read_csv(gwas_fn)
d = d[['SNP', 'Chromosome', 'Position ', 'P.value']]
d.columns = ['rs', 'chr', 'ps', 'p_wald']
# d = pd.read_csv(gwas_fn, sep='\t')
if bk.empty:
bk = d.copy()
bk.loc[bk.p_wald <= 1e-5, 'p_wald'] = 1e-5
for index, row in qtl.loc[qtl.phe_name == i, :].iterrows():
peak_pos = int(row['SNP'].split('_')[-1])
chrom = row['CHR']
sig_tmp = pd.concat([bk.loc[(bk.chr.astype(str) == str(chrom)) & (bk.ps >= peak_pos-1000000) & (bk.ps <= peak_pos+1000000), 'p_wald'],
d.loc[(bk.chr.astype(str) == str(chrom)) & (bk.ps >= peak_pos-1000000) & (bk.ps <= peak_pos+1000000), 'p_wald']], axis=1)
sig_tmp.columns = ['bk', 'phe']
bk.loc[(bk.chr.astype(str) == str(chrom)) & (bk.ps >= peak_pos-1000000) & (bk.ps <= peak_pos+1000000), 'p_wald'] = sig_tmp.apply(
lambda x: x['bk'] if x['bk'] < x['phe'] else x['phe'], axis=1)
bk = bk.loc[bk.p_wald <= 1e-2, :]
bk.loc[bk.p_wald <= 1e-20, 'p_wald'] = 1e-20
bk = bk[['rs', 'chr', 'ps', 'p_wald']]
thresholdi = robjects.FloatVector([1.0 / d.shape[0], 1e-6, 1e-5])
lim = -np.log10(min(bk['p_wald'])) + 2
bk.columns = ['SNP', 'Chromosome', 'Position', prefix]
base.sink('/dev/null')
robjects.r('source("'+utils_path+'/CMplot.r")')
#robjects.r('source("/home/debian/文档/MGWAP/compound_extract/plot/CMplot.r")')
CMplot = robjects.r['CMplot']
CMplot(bk, plot_type='m', col=robjects.StrVector(["grey30", "grey60"]), ylim=robjects.FloatVector([2, lim]),
threshold=thresholdi,
cex=robjects.FloatVector([0.5, 0.5, 0.5]), signal_cex=robjects.FloatVector([0.5, 0.5, 0.5]),
threshold_col=robjects.StrVector(['red', 'green', 'blue']), chr_den_col=robjects.rinterface.NULL,
amplify=True,
signal_pch=robjects.IntVector([19, 19, 19]), dpi=300,
signal_col=robjects.StrVector(['red', 'green', 'blue']), multracks=False, LOG10=True, file=t)
base.sink()
def qtl_anno(qtl, anno):
anno = anno[['geneid', 'position']]
anno.loc[:, 'chr'] = anno['position'].apply(lambda x: x.split(':')[0])
anno.loc[:, 'start'] = anno['position'].apply(lambda x: x.split(':')[1].split('-')[0])
anno.loc[:, 'end'] = anno['position'].apply(lambda x: x.split(':')[1].split('-')[1])
anno = anno[['chr', 'start', 'end', 'geneid']]
anno.columns = ['Chromosome', 'Start', 'End', 'geneid']
qtl_range = qtl[['CHR', 'qtl_start', 'qtl_end', 'phe_name']]
qtl_range.columns = ['Chromosome', 'Start', 'End', 'phe_name']
qtl_range = pr.PyRanges(qtl_range)
anno = pr.PyRanges(anno)
qtl_anno_intersect = qtl_range.join(anno, how='left')
qtl_anno = pd.DataFrame()
for k in sorted(qtl_anno_intersect.dfs.keys()):
qtl_anno = pd.concat([qtl_anno, qtl_anno_intersect.dfs[k]])
qtl_anno = qtl_anno[['Chromosome', 'Start', 'End', 'phe_name', 'geneid']]
qtl_anno.loc[:, 'Chromosome'] = qtl_anno.Chromosome.astype(str)
qtl_anno.loc[:, 'Start'] = qtl_anno.Start.astype(int)
qtl_anno.loc[:, 'End'] = qtl_anno.End.astype(int)
qtl.loc[:, 'CHR'] = qtl.CHR.astype(str)
qtl = pd.merge(qtl, qtl_anno, left_on=['CHR', 'qtl_start', 'qtl_end', 'phe_name'],
right_on=['Chromosome', 'Start', 'End', 'phe_name'])
qtl = qtl.drop(['Chromosome', 'Start', 'End'], axis=1)
qtl = qtl.groupby(['CHR', 'qtl_start', 'qtl_end', 'SNP', 'P', 'SP2_num', 'qtl_length', 'phe_name'])['geneid'].apply(';'.join).reset_index()
qtl.columns = ['CHR', 'qtl_start', 'qtl_end', 'SNP', 'P', 'SP2_num', 'qtl_length', 'phe_name', 'qtl_all_gene']
qtl.loc[qtl.qtl_all_gene == '-1', 'qtl_all_gene'] = 'nan'
return qtl
'''
Generate html report for SingleTrait
Author: CrazyHsu @ crazyhsu9527@gmail.com
Created on: 2020-08-26 20:33:52
Last modified: 2021-03-24 17:20:33
'''
################# Classes #################
class AllQtlStatistics():
def __init__(self):
self.totalSNPs = 0
self.qtlDetected = 0
self.medianQtlLen = 0
self.longestQtlLen = 0
self.shortestQtlLen = 0
self.totalAnnoGenes = 0
# self.totalTargetGenes = 0
# self.totalEnrichTargetGenes = 0
# self.aveEnrichGenes = 0
def getQtlLengthInfo(self, qtlLenSeries):
self.medianQtlLen = np.median(qtlLenSeries)
self.longestQtlLen = np.max(qtlLenSeries)
self.shortestQtlLen = np.min(qtlLenSeries)
def getTotalAnnoGenes(self, qtlAllGenesSeries):
self.totalAnnoGenes = len(self.mergeSeries2List(qtlAllGenesSeries))
# def getTotalTargetGenes(self, qtlAllMetaGenesSeries):
# self.totalTargetGenes = len(self.mergeSeries2List(qtlAllMetaGenesSeries))
# def getEnrichTargetGenes(self, myDataFrame):
# filteredSeries = myDataFrame.loc[myDataFrame.qtl_p_value <= 0.05, "qtl_meta_gene"]
# self.totalEnrichTargetGenes = len(self.mergeSeries2List(filteredSeries))
def getInit(self, tableData):
self.qtlDetected = len(tableData)
self.getQtlLengthInfo(tableData.qtl_length)
self.getTotalAnnoGenes(tableData.qtl_all_gene)
# self.getTotalTargetGenes(tableData.qtl_meta_gene)
# self.getEnrichTargetGenes(tableData)
def mergeSeries2List(self, mySeries, sep=";"):
return set().union(*mySeries.apply(lambda x: str(x).split(sep)).to_list())
class AllTraitStatistics():
def __init__(self):
self.totalTraits = 0
self.filteredTraits = 0
self.clusteredTraits = 0
self.unclusteredTraits = 0
def getFilteredTraits(self, traitSeries):
self.filteredTraits = len(self.mergeSeries2List(traitSeries, sep=","))
self.totalTraits = len(self.mergeSeries2List(traitSeries, sep=","))
def getClusteredTraits(self, labelSeries):
counter = Counter(labelSeries.to_list())
clustered = [i for i in counter if counter[i] != 1]
unclustered = [i for i in counter if counter[i] == 1]
self.clusteredTraits = len(clustered)
self.unclusteredTraits = len(unclustered)
# def getUnclusterTraits(self, labelSeries):
# self.clusteredTraits = len(set(labelSeries.to_list()))
def getInit(self, tableData):
self.getFilteredTraits(tableData.phe_name)
self.getClusteredTraits(tableData.phe_name)
# self.getUnclusterTraits(tableData.phe_name)
def mergeSeries2List(self, mySeries, sep=";"):
return set().union(*mySeries.apply(lambda x: str(x).split(sep)).to_list())
class SingleQtlStatistics():
def __init__(self, myRow, rowIndex):
self.qtlName = "QTL{}_chr{}-{}-{}".format(rowIndex, myRow.CHR, myRow.qtl_start, myRow.qtl_end)
self.qtlPosition = "{}:{}-{}".format(myRow.CHR, myRow.qtl_start, myRow.qtl_end)
self.peakSNP = "{}, {}".format(myRow.SNP, myRow.SNP)
self.traitNames = myRow.phe_name.split(",")
self.totalGenesInQtl = str(myRow.qtl_all_gene).split(";")
# self.targetGenesInQtl = str(myRow.qtl_meta_gene).split(";")
# self.enrichTargetGenes = 0
self.pvalue = myRow.P
################# Functions ##################
def resolveDir(dirName):
if not os.path.exists(dirName):
os.makedirs(dirName)
os.chdir(dirName)
def getGeneFunc(myGeneFuncFile, sep=","):
gene2func = pd.read_csv(myGeneFuncFile, sep=sep)
tmpDict = gene2func.to_dict("index")
geneFuncDict = {}
for i in tmpDict:
geneFuncDict[tmpDict[i]["geneid"]] = tmpDict[i]
return geneFuncDict
def getAllQtlSummary(allQtlStatistics, doc=None, tag=None, text=None, line=None):
with tag("div", id="qtlSummary"):
line("h1", "Summary of QTLs detected by local GWAS", style="text-align: center;")
multi_trait_File = "./Manhattan.multi_trait.jpg"
with tag("div", style="text-align: center;margin-top: 50px;"):
doc.stag("img", klass="img-fluid", src=multi_trait_File, alt="multi_trait")
with tag("div", style="margin-top: 100px;"):
with tag("div"):
line("h2", "QTL summarization criteria")
line("p", "The cutoffs in generating and filtering QTLs")
with tag("div"):
line("h2", "QTL statistics table")
with tag("div", klass="table-responsive"):
with tag("table", klass="table"):
with tag("thead"):
with tag("tr", klass="table-success"):
line("th", "Categories", style="width: 50%")
line("th", "Statistics value", style="width: 50%")
with tag("tbody"):
with tag("tr"):
line("td", "Genotype file")
line("td", "XXXXXX")
with tag("tr"):
line("td", "Number of SNPs for local GWAS")
line("td", str(allQtlStatistics.totalSNPs))
with tag("tr"):
line("td", "Number of detected QTLs")
line("td", str(allQtlStatistics.qtlDetected))
with tag("tr"):
line("td", "Median QTL length")
line("td", str(allQtlStatistics.medianQtlLen))
with tag("tr"):
line("td", "Longest QTL length")
line("td", str(allQtlStatistics.longestQtlLen))
with tag("tr"):
line("td", "Shortest QTL length")
line("td", str(allQtlStatistics.shortestQtlLen))
with tag("tr"):
line("td", "Total genes in the QTLs")
line("td", str(allQtlStatistics.totalAnnoGenes))
# with tag("tr"):
# line("td", "Total target genes in the QTLs")
# line("td", str(allQtlStatistics.totalTargetGenes))
# with tag("tr"):
# line("td", "Total enrichment of target genes")
# line("td", str(allQtlStatistics.totalEnrichTargetGenes))
# with tag("tr"):
# line("td", "Average enrichment of target genes")
# line("td", str(allQtlStatistics.aveEnrichGenes))
def getAllTraitSummary(allTraitStatistics, doc=None, tag=None, text=None, line=None):
# doc, tag, text, line = Doc().ttl()
with tag("div", id="traitSummary"):
line("h1", "Summary of omics traits detected by local GWAS", style="text-align: center;")
# heatmapFile = "/home/xufeng/xufeng/Projects/MODAS/xufeng1/assets/img/heatmap.png"
# with tag("div", style="text-align: center;margin-top: 50px;"):
# doc.stag("img", klass="img-fluid", src=heatmapFile, alt="heatmap")
with tag("div", style="margin-top: 100px;"):
with tag("div"):
line("h2", "Omics trait filtration criteria")
line("p", "Parameters and pipelines in filtering traits")
with tag("div"):
line("h2", "Trait statistics table")
with tag("div", klass="table-responsive"):
with tag("table", klass="table"):
with tag("thead"):
with tag("tr", klass="table-success"):
line("th", "Categories", style="width: 50%")
line("th", "Statistics value", style="width: 50%")
with tag("tbody"):
with tag("tr"):
line("td", "Omics trait type")
line("td", "Metabolome")
with tag("tr"):
line("td", "Total number of traits")
line("td", str(allTraitStatistics.totalTraits))
with tag("tr"):
line("td", "Filtered number of traits")
line("td", str(allTraitStatistics.filteredTraits))
with tag("tr"):
line("td", "Number of clustered traits")
line("td", str(allTraitStatistics.clusteredTraits))
# with tag("tr"):
# line("td", "Number of modules of clustered traits")
# line("td", "Cell 2")
with tag("tr"):
line("td", "Number of unclustered traits")
line("td", str(allTraitStatistics.unclusteredTraits))
# return doc.getvalue()
def getSingleQtlInfo(singleQtl, index, geneFuncDict, doc=None, tag=None, text=None, line=None):
with tag("div", klass="container"):
with tag("div", klass="row"):
with tag("div", klass="col"):
line("h1", "Summary in " + singleQtl.qtlName, style="text-align: center;")
with tag("div", style="margin-top: 100px;"):
with tag("div"):
line("h2", "QTL summarization criteria for whole -genome GWAS")
line("p", "Parameters and pipeline for filtering traits")
with tag("div"):
line("h2", "QTL statistics table")
with tag("div", klass="table-responsive"):
with tag("table", klass="table"):
with tag("thead"):
with tag("tr", klass="table-success"):
line("th", "QTL", style="width: 50%")
line("th", "Statistics", style="width: 50%")
with tag("tbody"):
with tag("tr"):
line("td", "QTL position")
line("td", str(singleQtl.qtlPosition))
with tag("tr"):
line("td", "Peak SNP ID and position")
line("td", str(singleQtl.peakSNP))
with tag("tr"):
line("td", "Number of total genes in the QTL")
line("td", str(len(singleQtl.totalGenesInQtl)))
# with tag("tr"):
# line("td", "Number of target genes in the QTL")
# line("td", str(len(singleQtl.targetGenesInQtl)))
# with tag("tr"):
# line("td", "Enrichment of target genes")
# line("td", str(singleQtl.enrichTargetGenes))
with tag("tr"):
line("td", "Enrichment significance vs background")
line("td", str(singleQtl.pvalue))
# with tag("div"):
# line("h2", "List of target genes in the QTL")
# with tag("div", klass="table-responsive"):
# with tag("table", klass="table"):
# with tag("thead"):
# with tag("tr", klass="table-success"):
# line("th", "Gene ID", style="width: 25%")
# line("th", "Alias ID", style="width: 25%")
# line("th", "Position", style="width: 25%")
# line("th", "Function", style="width: 25%")
# with tag("tbody"):
# for qtl in singleQtl.targetGenesInQtl:
# if qtl == "nan":
# continue
# with tag("tr"):
# line("td", str(geneFuncDict[qtl]["geneId"].strip()))
# line("td", str(geneFuncDict[qtl]["aliasId"].strip()))
# line("td", str(geneFuncDict[qtl]["position"].strip()))
# line("td", str(geneFuncDict[qtl]["function"].strip()))
with tag("div"):
line("h2", "List of total genes in the QTL")
with tag("div", klass="table-responsive"):
with tag("table", klass="table"):
with tag("thead"):
with tag("tr", klass="table-success"):
line("th", "Gene ID", style="width: 25%")
line("th", "Alias ID", style="width: 25%")
line("th", "Position", style="width: 25%")
line("th", "Function", style="width: 25%")
with tag("tbody"):
for qtl in singleQtl.totalGenesInQtl:
if qtl == "nan":
continue
with tag("tr"):
line("td", str(geneFuncDict[qtl]["geneid"].strip()))
line("td", str(geneFuncDict[qtl]["aliasid"].strip()))
line("td", str(geneFuncDict[qtl]["position"].strip()))
line("td", str(geneFuncDict[qtl]["function"].strip()))
# def getSingleTrait(traitName, doc=None, tag=None, text=None, line=None):
# with tag("div", klass="container"):
# with tag("div", klass="row"):
# with tag("div", klass="col"):
# line("h1", "Details in " + traitName, style="text-align: center;")
# doc.stag("img", klass="img-fluid", src="../../assets/img/manhattan.jpg")
# with tag("div", klass="row"):
# with tag("div", klass="col-6"):
# doc.stag("img", klass="img-fluid", src="../../assets/img/qqplot.png")
# with tag("div", klass="col-6"):
# doc.stag("img", klass="img-fluid", src="../../assets/img/boxplot.png")
def getListItem(data, qtlName=None, traitName=None, doc=None, tag=None, text=None, line=None, mainPage=False):
for index, row in data.iterrows():
qtlItem = SingleQtlStatistics(row, index)
if qtlName and qtlName == qtlItem.qtlName:
expand = "true"
faPlusOrMinus = "fa-minus"
myClass = "list-unstyled collapse nav nav-pills show"
active = " active"
else:
expand = "false"
faPlusOrMinus = "fa-plus"
myClass = "list-unstyled collapse nav nav-pills"
active = ""
if mainPage:
relativeDir = os.path.join("", qtlItem.qtlName)
else:
if qtlName == qtlItem.qtlName:
relativeDir = ""
else:
relativeDir = os.path.join("../", qtlItem.qtlName)
with tag("li"):
with tag("div", klass="qtlItem" + active):
with tag("a", ("href", os.path.join(relativeDir, qtlItem.qtlName + ".html")), klass="qtlLink"):
text(qtlItem.qtlName)
with tag("a", ("href", "#" + qtlItem.qtlName), ("data-toggle", "collapse"), ("aria-expanded", expand)):
line("i", "", klass="fa " + faPlusOrMinus)
with tag("ul", ("class", myClass), ("id", qtlItem.qtlName), ("aria-expanded", expand)):
for i in qtlItem.traitNames:
with tag("li"):
href = os.path.join(relativeDir, i + ".html")
if traitName and traitName == i:
with tag("a", ("href", href), ("class", "active"), ("aria-selected", "true")):
line("i", "", klass="fa fa-link")
text(" " + i)
else:
with tag("a", ("href", href), ("aria-selected", "false")):
line("i", "", klass="fa fa-link")
text(" " + i)
def generateMainPage(data, allQtlStatistics, allTraitStatistics):
doc, tag, text, line = Doc().ttl()
doc.asis('<!DOCTYPE html>')
with tag('html'):
with tag('head'):
doc.stag('meta', charset='utf-8')
doc.stag('meta', name='viewport', content='width=device-width, initial-scale=1.0, shrink-to-fit=no')
line('title', 'MODAS main page')
doc.stag('link', rel='stylesheet', href='assets/bootstrap/css/bootstrap.min.css')
doc.stag('link', rel='stylesheet', href='assets/fonts/font-awesome.min.css')
doc.stag('link', rel='stylesheet', href="assets/css/modas.css")
# doc.stag('link', rel='stylesheet', href="assets/css/styles.css")
with tag("body"):
with tag("div", id="sidebar-test"):
with tag("div", klass="sidebar-header"):
with tag("h2"):
line("a", "MODAS", href="mainPage.html", klass="modas")
with tag("ul"):
getListItem(data, doc=doc, tag=tag, text=text, line=line, mainPage=True)
with tag("div", klass="content"):
with tag("div", klass="container"):
with tag("div", klass="row"):
with tag("div", klass="col"):
getAllQtlSummary(allQtlStatistics, doc=doc, tag=tag, text=text, line=line)
doc.stag("hr", style="margin-bottom: 50px;margin-top: 50px;")
getAllTraitSummary(allTraitStatistics, doc=doc, tag=tag, text=text, line=line)
line("script", "", src="assets/js/jquery.min.js")
line("script", "", src="assets/bootstrap/js/bootstrap.min.js")
line("script", "", src="assets/js/modas.js")
mainPageOut = open("mainPage.html", "w")
res = indent(doc.getvalue(), indentation=" ")
print(res, file=mainPageOut)
mainPageOut.close()
def generateSingleQtlPage(data, geneFuncDict):
for index, row in data.iterrows():
qtlItem = SingleQtlStatistics(row, index)
if not os.path.exists(qtlItem.qtlName):
os.makedirs(qtlItem.qtlName)
out = open(os.path.join(qtlItem.qtlName, qtlItem.qtlName + ".html"), "w")
doc, tag, text, line = Doc().ttl()
doc.asis('<!DOCTYPE html>')
with tag('html'):
with tag('head'):
doc.stag('meta', charset='utf-8')
doc.stag('meta', name='viewport', content='width=device-width, initial-scale=1.0, shrink-to-fit=no')
line('title', 'Summary information in QTL ' + qtlItem.qtlName)
doc.stag('link', rel='stylesheet', href='../assets/bootstrap/css/bootstrap.min.css')
doc.stag('link', rel='stylesheet', href='../assets/fonts/font-awesome.min.css')
doc.stag('link', rel='stylesheet', href="../assets/css/modas.css")
# doc.stag('link', rel='stylesheet', href="assets/css/styles.css")
with tag("body"):
with tag("div", id="sidebar-test"):
with tag("div", klass="sidebar-header"):
with tag("h2"):
line("a", "MODAS", href="../mainPage.html", klass="modas")
with tag("ul"):
getListItem(data, qtlName=qtlItem.qtlName, doc=doc, tag=tag, text=text, line=line)
with tag("div", klass="content"):
getSingleQtlInfo(qtlItem, index, geneFuncDict, doc=doc, tag=tag, text=text, line=line)
line("script", "", src="../assets/js/jquery.min.js")
line("script", "", src="../assets/bootstrap/js/bootstrap.min.js")
line("script", "", src="../assets/js/modas.js")
customJs = '''
<script>
var offestFromTop = %d * 45 + 68;
$('#sidebar-test').scrollTop(offestFromTop);
function clickItem(event) {
var target = event.currentTarget;
$(target).parent().removeClass(".active").addClass(".active");
var index = $("div.qtlItem").index($(this).parent());
var offestFromTop = index * 45 + 68;
$('#sidebar-test').scrollTop(offestFromTop);
}
if ($("div.qtlItem .qtlLink")) {
var qtlLink = $("div.qtlItem .qtlLink");
for (var i = 0; i < qtlLink.length; i++) {
var item = qtlLink[i];
item.onclick = clickItem;
}
}
</script>
''' % (index)
doc.asis(customJs)
res = indent(doc.getvalue(), indentation=" ")
print(res, file=out)
out.close()
def generateSingleTraitPage(data, manhattanDir, qqDir, boxplotDir):
for index, row in data.iterrows():
qtlItem = SingleQtlStatistics(row, index)
for traitName in qtlItem.traitNames:
# getListItem(qtlItem)
out = open(os.path.join(qtlItem.qtlName, traitName + ".html"), "w")
doc, tag, text, line = Doc().ttl()
doc.asis('<!DOCTYPE html>')
with tag('html'):
with tag('head'):
doc.stag('meta', charset='utf-8')
doc.stag('meta', name='viewport', content='width=device-width, initial-scale=1.0, shrink-to-fit=no')
line('title', 'Detailed information in trait ' + traitName)
doc.stag('link', rel='stylesheet', href='../assets/bootstrap/css/bootstrap.min.css')
doc.stag('link', rel='stylesheet', href='../assets/fonts/font-awesome.min.css')
doc.stag('link', rel='stylesheet', href="../assets/css/modas.css")
# doc.stag('link', rel='stylesheet', href="assets/css/styles.css")
with tag("body"):
with tag("div", id="sidebar-test"):
with tag("div", klass="sidebar-header"):
with tag("h2"):
line("a", "MODAS", href="../mainPage.html", klass="modas")
with tag("ul"):
getListItem(data, qtlItem.qtlName, traitName, doc=doc, tag=tag, text=text, line=line)
manhattanFile = glob.glob(os.path.join(manhattanDir, "Manhattan.*_{}.jpg".format(traitName)))[0]
boxplotFile = glob.glob(os.path.join(boxplotDir, "{}_*.jpg".format(traitName)))[0]
qqplotFile = glob.glob(os.path.join(qqDir, "QQplot.*_{}.jpg".format(traitName)))[0]
with tag("div", klass="content"):
with tag("div", klass="container"):
with tag("div", klass="row"):
with tag("div", klass="col"):
line("h1", "Details in " + traitName, style="text-align: center;")
doc.stag("img", klass="img-fluid", src='../'+'/'.join(manhattanFile.split('/')[-2:]))
with tag("div", klass="row"):
with tag("div", klass="col-6"):
doc.stag("img", klass="img-fluid", src='../'+'/'.join(qqplotFile.split('/')[-2:]))
with tag("div", klass="col-6"):
doc.stag("img", klass="img-fluid", src='../'+'/'.join(boxplotFile.split('/')[-2:]))
line("script", "", src="../assets/js/jquery.min.js")
line("script", "", src="../assets/bootstrap/js/bootstrap.min.js")
line("script", "", src="../assets/js/modas.js")
customJs = '''
<script>
var offestFromTop = %d * 45 + 68;
$('#sidebar-test').scrollTop(offestFromTop);
function clickItem(event) {
var target = event.currentTarget;
$(target).parent().removeClass(".active").addClass(".active");
var index = $("div.qtlItem").index($(this).parent());
var offestFromTop = index * 45 + 68;
$('#sidebar-test').scrollTop(offestFromTop);
}
if ($("div.qtlItem .qtlLink")) {
var qtlLink = $("div.qtlItem .qtlLink");
for (var i = 0; i < qtlLink.length; i++) {
var item = qtlLink[i];
item.onclick = clickItem;
}
}
</script>
''' % (index)
doc.asis(customJs)
res = indent(doc.getvalue(), indentation=" ")
print(res, file=out)
out.close()
def generateHtml(qtl_anno, myFuncFile, out_dir, totalSNPs):
manhattanDir = os.path.abspath(out_dir+'/manhattan_plot')
qqplotDir = os.path.abspath(out_dir+'/qqplot')
boxplotDir = os.path.abspath(out_dir+'/boxplot')
allQtlStatistics = AllQtlStatistics()
allQtlStatistics.getInit(qtl_anno)
allQtlStatistics.totalSNPs = totalSNPs
allTraitStatistics = AllTraitStatistics()
allTraitStatistics.getInit(qtl_anno)
geneFuncDict = getGeneFunc(myFuncFile, "\t")
resolveDir(out_dir)
generateMainPage(qtl_anno, allQtlStatistics, allTraitStatistics)
generateSingleQtlPage(qtl_anno, geneFuncDict)
generateSingleTraitPage(qtl_anno, manhattanDir, qqplotDir, boxplotDir)
Datasets
Models
