functions{ #include /functions.stan vector f(vector alpha, vector beta , real[] xr , int[] xi ){ int M = size(xr); int N_DNA = xi[M+1, 1]; int N_RNA = xi[M+1, 2]; int a[M,N_DNA] = xi[1:M,1:N_DNA]; int b[M,(N_DNA+1):(2*N_DNA)] = xi[1:M,(N_DNA+1):(2*N_DNA)]; int k[M, M,(2*N_DNA+1):(2*N_DNA+N_RNA)] = xi[1:M,(2*N_DNA+1):(2*N_DNA+N_RNA)]; int m[M, (2*N_DNA+N_RNA+1):2*(N_DNA+N_RNA)] = xi[1:M,(2*N_DNA+N_RNA+1):2*(N_DNA+N_RNA)]; int p[M] = alpha[1,1:M]; int theta[M] = alpha[2,1:M]; int q[M] = alpha[3,1:M]; int qi[M,(N_RNA+2)-3] = alpha[4:(N_RNA+2), 1:M]'; for(j in 1:N_VARIANTS) { p[j] ~ betaModeConc(v[j],beta[2]); log(theta[j])/beta[3] ~ normal(0,1); // theta~lognormal(0,s) l +=-log(theta[j])-log(beta[3]); // Jacobian for(i in 1:N_RNA) qi[j,i] ~ betaModeConc(q[j],beta[1]); for(i in 1:N_DNA) a[j,i] ~ binomial(a[j,i]+b[j,i],p[j]); for(i in 1:N_RNA) k[j,i] ~ binomial(k[j,i]+m[j,i],qi[j,i]); } return [l]; } } data { int N_VARIANTS; // number of variants in the model real v[N_VARIANTS]; // alt allele freq in VCF file int N_DNA; // number of DNA replicates int a[N_VARIANTS,N_DNA]; // DNA alt read counts int b[N_VARIANTS,N_DNA]; // DNA ref read counts int N_RNA; // number of RNA replicates int k[N_VARIANTS,N_RNA]; // RNA alt read counts int m[N_VARIANTS,N_RNA]; // RNA ref read counts } transformed data { int N_THREAD = 2; int M = N_VARIANTS/N_THREAD; vector[2] Number; //Define a vector to include N_DNA, N_RNA Numebr[1] = N_DNA; //include N_DNA Number[2] = N_RNA; //include N_RNA int xi[M+1, 2*(N_DNA + N_RNA), N_THREAD]; //extra 1 row is to store Number real xr[N_THREAD, M]; for ( i in 1:N_THREAD ) { int j = 1 + (i-1)*M; int k = i*M; xi[1:M, 1:N_DNA, i] = a[j:k,N_DNA]; xi[1:M, N_DNA+1:2*N_DNA, i] = b[j:k,N_DNA]; xi[1:M, 2*N_DNA+1:2*N_DNA+N_RNA, i] = k[j:k,N_RNA]; xi[1:M,2*N_DNA+N_RNA+1:2*(N_DNA + N_RNA), i] = m[j:k,N_RNA]; xr[i] = v[j:k]; xi[M+1, 1:2, i] = Number[1:2]; //store the vector include N_DNA, N_RNA } } parameters { real p[N_VARIANTS]; // alt allele freq in DNA real qi[N_VARIANTS,N_RNA]; // alt freqs in RNA real theta[N_VARIANTS]; real c1; // concentration parameter of beta prior for qi real c2; // concentration parameter of beta prior for p real s; // std.dev. parameter of lognormal prior for theta vector[4] beta[N_THREAD]; ////Universal Parameter for 2 shards for (i in 1:N_THREAD){ beta[i, 1] = c1; beta[i, 2] = c2; beta[i, 3] = s; } matrix[N_RNA+3, M] alpha[N_THREAD] for (i in 1:N_THREAD){ int j = 1 + (i-1)*M; int k = i*M; alpha[i, 1, 1:M] = p[j:k]; alpha[i, 2, 1:M] = theta[j:k]; alpha[i, 4:(N_RNA+2), 1:M] = qi[j:k,1:N_RNA]'; } } transformed parameters { // ORDER MATTERS! real q[M]; // alt allele freq in RNA for(j in 1:M) q[j]=theta[j]*p[j]/(1.0-p[j]+theta[j]*p[j]); ///////////////////////////////// for (i in 1:N_THREAD){ int j = 1 + (i-1)*M; int k = i*M; alpha[i, 3, 1:M] = q[j:k]; } } model { // Parameters c1 ~ gamma(1.1, 0.005); c2 ~ gamma(1.1, 0.005); s ~ gammaModeSD(1,1); target += sum(map_rect(f, beta , alpha , xr , xi)); }