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