Hi I have gene expression, phenotypic (sex, age etc info) and snps information in one file (RA_122) and gene snp pair in another (df). I trying to run a linear mixed effects model using python to generate r codes as you can see below FOR THE EQTL ANALYSIS.

def R(cmd):
    print(cmd)

    import pandas

    R('library("nlme")')
    R('RA_122=read.table("Gene_phenotype_snps.txt",sep="t")')
    R('sink("Test.txt")')

    df = pandas.read_csv('SNPs_gene.txt', sep='t')

    for snp, gene in zip(df['Gene_stable_id'], df['snp']):
        R('formula= %s ~ %s * SE.Allele + Age + RIN + Sex'%(snp, gene))
        R(r'''tryCatch({test.lme.ptpn22=lme(formula, RA_122, random = ~ 1 | IDexpressionset, na.action=na.omit)
        effects=summary(test.lme.ptpn22)$tTable[,"Value"]
        pvalues=summary(test.lme.ptpn22)$tTable[,"p-value"]
        cat(c("%s", "%s", pvalues, effects)); cat("n")}, error=function(e){})'''%(gene, snp))
        R('sink()')

When I run this code it as

eQtl_RA.py > R_codes.txt

It successfully generates a list of R codes for hundred of snp/genes that look like this in the file R_codes.txt.

library("nlme")
ACPA.122snps=read.table("Final_CD4_factor.txt",sep="    ")
sink("Test-eQTLs_CD4.txt")
formula= ENSG00000158828 ~ rs3121680 * SE.Allele + Age + RIN + Sex
tryCatch({test.lme.ptpn22=lme(formula, ACPA.122snps, random = ~ 1 | IDexpressionset, na.action=na.omit)
effects=summary(test.lme.ptpn22)$tTable[,"Value"]
pvalues=summary(test.lme.ptpn22)$tTable[,"p-value"]
cat(c("rs3121680", "ENSG00000158828", pvalues, effects)); cat("n")}, error=function(e){})
formula= ENSG00000132849 ~ rs4456122 * SE.Allele + Age + RIN + Sex
tryCatch({test.lme.ptpn22=lme(formula, ACPA.122snps, random = ~ 1 | IDexpressionset, na.action=na.omit)
effects=summary(test.lme.ptpn22)$tTable[,"Value"]
pvalues=summary(test.lme.ptpn22)$tTable[,"p-value"]
cat(c("rs4456122", "ENSG00000132849", pvalues, effects)); cat("n")}, error=function(e){})
formula= ENSG00000081026 ~ rs10858000 * SE.Allele + Age + RIN + Sex
tryCatch({test.lme.ptpn22=lme(formula, ACPA.122snps, random = ~ 1 | IDexpressionset, na.action=na.omit)
effects=summary(test.lme.ptpn22)$tTable[,"Value"]
pvalues=summary(test.lme.ptpn22)$tTable[,"p-value"]
cat(c("rs10858000", "ENSG00000081026", pvalues, effects)); cat("n")}, error=function(e){})
formula= ENSG00000081026 ~ rs11102648 * SE.Allele + Age + RIN + Sex
tryCatch({test.lme.ptpn22=lme(formula, ACPA.122snps, random = ~ 1 | IDexpressionset, na.action=na.omit)
effects=summary(test.lme.ptpn22)$tTable[,"Value"]
pvalues=summary(test.lme.ptpn22)$tTable[,"p-value"]
cat(c("rs11102648", "ENSG00000081026", pvalues, effects)); cat("n")}, error=function(e){})
.
.
.

But the results do not sink into the Test.txt file automatically. And when I run

 Rscript  R_codes.txt

I get

Null
Null
Null 

for hundred of lines in the Test.txt file instead of results.

Any suggesstions on what could be wrong? My data has been called as.factors when I prepared that files earlier.

Regards



Source link