In [1]:
%cd ../data
In [2]:
import numpy as np
In [3]:
def clip(s):
return max(0,int(float(s)+.5))
In [4]:
from csv import reader, excel_tab
orfs = []
names = []
data = []
fin = reader(open("GSE88801_kallisto_est_counts_thresh1.cdt"),dialect=excel_tab)
header = next(fin)[4:]
eweights = next(fin)[4:]
for row in fin:
orfs.append(row[1])
names.append(row[2])
data.append([clip(i) for i in row[4:]])
In [5]:
C = np.array(data)
C.shape
Out[5]:
In [6]:
from csv import reader
In [7]:
fp = reader(open("sample_table_v2.csv"))
sample_header = next(fp)
samples = []
name2sample = {}
for i in fp:
samples.append(i)
name2sample[i[0]] = i
In [8]:
sample_header
Out[8]:
In [9]:
out = open("sample_table_v2.tdt","w")
out.write("\t".join(sample_header+["state"])+"\n")
for i in header:
s = name2sample[i]
out.write("\t".join(s + ["%s.%s.%s" % (s[2],s[1],s[3])])+"\n")
out.close()
In [10]:
out = open("GSE88801_kallisto_est_counts_thresh1.txt","w")
out.write("\t".join(["gene"]+header)+"\n")
for (name,row) in zip(names,C):
out.write("\t".join([name]+[str(i) for i in row])+"\n")
out.close()
In [11]:
%load_ext rpy2.ipython
In [12]:
%%R
library(limma)
library(edgeR)
In [13]:
%%R
samples <- read.delim("sample_table_v2.tdt")
print(summary(samples))
In [14]:
%%R
state <- samples$state
d <- model.matrix(~0+state)
colnames(d) <- gsub("state","",colnames(d))
print(colnames(d))
In [15]:
%%R
C <- read.delim("GSE88801_kallisto_est_counts_thresh1.txt",row.names=1)
dge <- DGEList(counts=C)
dge <- calcNormFactors(dge)
v <- voom(dge, d, plot = TRUE)
In [16]:
%%R
fit <- lmFit(v, d)
contrast.matrix <- makeContrasts(# S1 contrast
BMDM.uninfected.24-BMDM.uninfected.4,
# S2 contrasts
BMDM.Live.4-BMDM.uninfected.4,BMDM.Dead.4-BMDM.uninfected.4,
BMDM.Live.24-BMDM.uninfected.24,BMDM.Dead.24-BMDM.uninfected.24,
J774.Live.4-J774.uninfected.4,J774.Dead.4-J774.uninfected.4,
J774.Live.24-J774.uninfected.24,J774.Dead.24-J774.uninfected.24,
levels=d)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
print(summary(decideTests(fit2)))
In [17]:
%%R
print(summary(decideTests(fit2,lfc=1)))
In [18]:
%%R
print(summary(decideTests(fit2,lfc=.27)))
In [19]:
%%R
print(topTable(fit2))
In [20]:
%%R
print(colnames(fit2$coefficients))
In [21]:
%%R
print(topTable(fit2, coef="BMDM.uninfected.24 - BMDM.uninfected.4", lfc=1, p.value = .05))
In [22]:
%%R
for(tc in colnames(fit2$coefficients)){
print(tc)
# Extract all genes significantly differential on this contrast for a 2x fold change cutoff and 5% FDR
# Use write.csv rather than write.table for clean compatibility with python's csv.reader
write.csv(topTable(fit2, coef=tc, n = 40000, lfc=1, p.value = .05),
paste("limma1.",gsub(" ","",tc),".t0.csv",sep=""))
# Extract the adjusted p-values for this contrast for all genes, independent of significance
write.csv(topTable(fit2, coef=tc, n = 40000),
paste("limma1.",gsub(" ","",tc),".t1.csv",sep=""))
}
In [ ]: