Skip to content

Commit ae24b7e

Browse files
Jihong528Jihong528
committed
First Commit
0 parents  commit ae24b7e

34 files changed

+10202
-0
lines changed

‎BetaBinom.R

+141
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
library('VGAM')
2+
# library('rjags')
3+
# library('overlapping')
4+
library('aod')
5+
options(scipen=999) # scipen 科学计数方法 小于999的都会显示出来
6+
7+
args = commandArgs(trailingOnly=TRUE)
8+
filename = args[1]
9+
count = read.table(paste(filename, '.count.txt',sep=''))
10+
N = as.numeric(colSums(count))
11+
tissue=c(rep('CSF',8),rep('PBMC',10))
12+
states = c('MS','MS','MS','MS','control','control','control','control',
13+
'MS','MS','MS','MS','MS','control','control','control','control','control')
14+
15+
16+
# BetaBinomial1 = "
17+
# model{
18+
# for(i in 1:N){
19+
# x[i] ~ dbinom(p1[i]*c[i]+p0[i]*(1-c[i]),n[i])
20+
# p0[i] ~ dbeta(a0,b0) T(,1-1e-100)
21+
# p1[i] ~ dbeta(a1,b1) T(,1-1e-100)
22+
# }
23+
# a0 ~ dgamma(1e-10, 1e-10)
24+
# b0 ~ dgamma(1e-10, 1e-10)
25+
# a1 ~ dgamma(1e-10, 1e-10)
26+
# b1 ~ dgamma(1e-10, 1e-10)
27+
# }"
28+
29+
res1 = BetaBinomialReg(count,N,tissue=='CSF',paste(filename,'.betabinomreg.CSF.csv',sep=''))
30+
31+
BetaBinomialReg = function(count,N,C,filename){
32+
res = lapply(c(1:length(count[,1])),function(i){
33+
x = as.numeric(count[i,])
34+
# fit = vglm(t(rbind( x, N-x )) ~ C, betabinomialff)
35+
# p = coef(summary(fit))[c(3,4),4]
36+
# pred = fitted(fit)
37+
# fc = pred[C==T][1] / pred[C==F][1]
38+
data = data.frame(x=x,N=N,C=as.factor(C))
39+
fm1 <- betabin(cbind(x, N - x) ~ C, ~ 1, data)
40+
res = summary(fm1)@Coef
41+
p = res[2,4]
42+
# pred = fitted(fm1)
43+
# fc = pred[C==T][1] / pred[C==F][1]
44+
fc = mean(x[C==T]/N[C==T])/mean(x[C==F]/N[C==F])
45+
return(c(p,fc))
46+
})
47+
res = do.call(cbind,res)
48+
write.table(t(res),col.names=F,row.names=F,quote=F,sep=',',file=filename)
49+
return(res)
50+
}
51+
52+
53+
54+
# BinomialReg = function(count,N,C,filename){
55+
# res = lapply(c(1:length(count[,1])),function(i){
56+
# x = as.numeric(count[i,])
57+
# df <- data.frame(
58+
# TotalCells = N,
59+
# numInCluster = x,
60+
# Genotype = as.factor(C)
61+
# )
62+
# result <- glm(numInCluster/TotalCells~Genotype,
63+
# data=df, family=binomial,
64+
# weights=TotalCells)
65+
# p = summary(result)$coef[2,4]
66+
# return(c(p,mean(x[C==T]/N[C==T])/mean(x[C==F]/N[C==F]) ))
67+
# })
68+
# res = do.call(cbind,res)
69+
# write.table(t(res),col.names=F,row.names=F,quote=F,sep=',',file=filename)
70+
# return(res)
71+
# }
72+
73+
# ComputeDP = function(count,N,C,filename){
74+
# pdf(paste(filename,'.pdf',sep=''))
75+
# OV = lapply(c(1:length(count[,1])),function(i){
76+
# X = count[i,]
77+
# model <- jags.model(textConnection(BetaBinomial1),
78+
# data = list('x' = as.numeric(X),
79+
# 'n' = as.numeric(N),
80+
# 'c' = C,
81+
# 'N' = length(N)))
82+
# update(model, 5000)
83+
# params <- jags.samples(model,c('a0','a1','b0','b1'),50000)
84+
# mu0 = params[[1]][1,,1]/(params[[1]][1,,1]+params[[3]][1,,1])
85+
# mu1 = params[[2]][1,,1]/(params[[2]][1,,1]+params[[4]][1,,1])
86+
# par(mfrow=c(3,3))
87+
# lower = min(c(mu0,mu1))
88+
# upper = max(c(mu0,mu1))
89+
# step = (upper-lower)/100
90+
# hist(mu0,breaks=seq(lower,upper,step),col='blue',border='blue')
91+
# hist(mu1,breaks=seq(lower,upper,step),col='red',add=T,border='red')
92+
# res = c(overlap(list(log(mu0),log(mu1)))$OV, mean(mu1)/mean(mu0))
93+
# return(list(res,mu0,mu1))
94+
# })
95+
# mu0 = lapply(OV,function(X){X[[2]]})
96+
# mu1 = lapply(OV,function(X){X[[3]]})
97+
# OV = sapply(OV,function(X){X[[1]]})
98+
# plot(y=-log(OV[1,]+1e-10,10),x=log(OV[2,],10))
99+
# dev.off()
100+
# write.table(t(OV),col.names=F,row.names=F,quote=F,sep=',',file=filename)
101+
# return(list(OV,mu0,mu1))
102+
# }
103+
104+
# OV1 = ComputeDP(count,N,tissue=='CSF','Bayesian.CSF.csv')
105+
# OV2 = ComputeDP(count,N,states=='MS','Bayesian.MS.csv')
106+
# OV3 = ComputeDP(count[,tissue=='CSF'],N[tissue=='CSF'],states[tissue=='CSF']=='MS','Bayesian.CSFMS.csv')
107+
# OV4 = ComputeDP(count[,tissue=='PBMC'],N[tissue=='PBMC'],states[tissue=='PBMC']=='MS','Bayesian.PBMCMS.csv')
108+
# OV5 = ComputeDP(count[,states=='control'],N[states=='control'],tissue[states=='control']=='CSF','Bayesian.CSF_control.csv')
109+
110+
# tres = c(4.29382293, 3.33025022, 7.61345436, 0.97755275, 0.91960306,
111+
# 1.5547467 , 0.36301796, 0.78215611, 0.57243582, 0.46252548,
112+
# 5.18488324, 2.8403952 , 0.78215611, 0.59855511, 2.02233809,
113+
# 0.91651382, 0.46077512, 1.2522841 , 2.47588862, 0.21768479,
114+
# 0.71656765, 8.59107891)
115+
116+
# count[12, states=='MS' & tissue=='CSF']
117+
# res1 = BinomialReg(count,N,tissue=='CSF','count_test/binomreg.CSF.csv')
118+
# res2 = BinomialReg(count,N,states=='MS','count_test/binomreg.MS.csv')
119+
# res3 = BinomialReg(count[,tissue=='CSF'],N[tissue=='CSF'],states[tissue=='CSF']=='MS','count_test/binomreg.CSFMS.csv')
120+
# res4 = BinomialReg(count[,tissue=='PBMC'],N[tissue=='PBMC'],states[tissue=='PBMC']=='MS','count_test/binomreg.PBMCMS.csv')
121+
# res5 = BinomialReg(count[,states=='control'],N[states=='control'],tissue[states=='control']=='CSF','count_test/binomreg.CSF_control.csv')
122+
123+
124+
res1 = BetaBinomialReg(count,N,tissue=='CSF',paste(filename,'.betabinomreg.CSF.csv',sep=''))
125+
res2 = BetaBinomialReg(count,N,states=='MS',paste(filename,'.betabinomreg.MS.csv',sep=''))
126+
res3 = BetaBinomialReg(count[,tissue=='CSF'],N[tissue=='CSF'],states[tissue=='CSF']=='MS',
127+
paste(filename,'.betabinomreg.CSFMS.csv',sep=''))
128+
res4 = BetaBinomialReg(count[,tissue=='PBMC'],N[tissue=='PBMC'],states[tissue=='PBMC']=='MS',
129+
paste(filename,'.betabinomreg.PBMCMS.csv',sep=''))
130+
res5 = BetaBinomialReg(count[,states=='control'],N[states=='control'],tissue[states=='control']=='CSF',
131+
paste(filename,'.betabinomreg.CSF_control.csv',sep=''))
132+
133+
134+
135+
# 晓竹师姐调试出来的代码
136+
count = read.csv('female_sum_neuron_celltype.csv',header = T, row.names = 1)
137+
countnew = t(count)
138+
N = as.numeric(colSums(countnew))
139+
states = c('GABA1','GABA2','GABA3','GABA4','GABA5','GABA6','Glu1','Glu2','Glu3')
140+
filename="neuron"
141+
res1gaba1 = BetaBinomialReg(countnew,N,states=='GABA1',paste(filename,'.betabinomreg.CSF.csv',sep=''))

‎BetaBinom.R.html

+1,622
Large diffs are not rendered by default.

‎CalculateReadCountBarcode.R

+16
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
CalculateReadCountBarcode <- function(mol.info.file, row.prefix){
2+
3+
library(Seurat)
4+
library(DropletUtils)
5+
molinfo <- read10xMolInfo(mol.info.file)
6+
selected_info <- molinfo$data[1:5]
7+
selected_info$cell <- factor(selected_info$cell)
8+
cell_read_sum <- tapply(selected_info$reads,selected_info$cell,sum)
9+
10+
cell_read_sum <- as.matrix(cell_read_sum)
11+
rownames(cell_read_sum) <- paste0(row.prefix, rownames(cell_read_sum))
12+
colnames(cell_read_sum) <- "reads_sum"
13+
print(dim(cell_read_sum))
14+
print(head(cell_read_sum))
15+
return(cell_read_sum)
16+
}

0 commit comments

Comments
 (0)