<jupyter_start><jupyter_text>mpra_analyze_SL<NAME>12/10/2020<jupyter_code>library(MPRAnalyze)
library(tidyverse)
library(pheatmap)
# library(mpra)
# library(QuASAR)
# library(BiocParallel)
# library(stringr)
# library(reshape2)
# library(ggplot2)
# library(gridExtra)
# library(VennDiagram)

getwd()
list.files(getwd())
count_rna = read.csv('dedup_counts/count_matrix.csv',row.names=1)
head(count_rna)
colnames(count_rna)
get_tables = function(miseq_out){
miseq_out$id = as.character(miseq_out$id)
miseq_out = miseq_out%>%
    select(id, count)%>%
    separate(id, c("chr", "loc", "allele_info", "letter", "bc_id"), sep='[.]')%>%
    mutate(allele = ifelse(grepl("Ref",miseq_out$id),'ref','alt'))%>%
    mutate(column = str_c(allele, bc_id, sep = '_'))%>%
    mutate(snp_name = str_c(chr, loc, sep = '_'))%>%
    select(snp_name, count, column, allele, bc_id )%>%
    group_by(snp_name, column)%>%
    summarise(count = sum(count),
            allele = unique(allele),
            bc_id = unique(bc_id))%>%
    ungroup()

    colAnnon = miseq_out%>%
        select( -snp_name,-count)%>%
        mutate(barcode_allelic =  as.factor(str_c(allele, bc_id, sep = '_')))%>%
        distinct()%>%
        arrange(column)%>%
        column_to_rownames('column')
    
    countTable = miseq_out%>%
        select(-allele, -bc_id)%>%
        spread(column,count)%>%
        replace(is.na(.), 0)%>%
        arrange(snp_name)%>%
        column_to_rownames('snp_name')
#     countTable = as.matrix(countTable)
    return(list(colAnnon, countTable))
}

get_tables_withbatch = function(miseq_out){
miseq_out$id = as.character(miseq_out$id)
miseq_out = miseq_out%>%
    select(id, count,batch)%>%
    separate(id, c("chr", "loc", "allele_info", "letter", "bc_id"), sep='[.]')%>%
    mutate(allele = ifelse(grepl("Ref",miseq_out$id),'ref','alt'))%>%
    mutate(column = str_c(allele, bc_id, batch, sep = '_'))%>%
    mutate(snp_name = str_c(chr, loc, sep = '_'))%>%
    select(snp_name, count, column, allele, bc_id, batch )%>%
    group_by(snp_name, column, batch)%>%
    summarise(count = sum(count),
            allele = unique(allele),
            bc_id = unique(bc_id))%>%
    ungroup()

    colAnnon = miseq_out%>%
        select( -snp_name,-count)%>%
        mutate(barcode_allelic =  as.factor(str_c(allele, bc_id, sep = '_')))%>%
        distinct()%>%
        arrange(column)%>%
        column_to_rownames('column')

    countTable = miseq_out%>%
        select(-allele, -bc_id, -batch)%>%
        spread(column,count)%>%
        replace(is.na(.), 0)%>%
        arrange(snp_name)%>%
        column_to_rownames('snp_name')
#     countTable = as.matrix(countTable)

    return(list(colAnnon, countTable))
}

# read_rna_count_csv = function(rna_count_csv){
    
# }
# read in dna
dna_df = read.csv('../miseq_mg_081019/libpsych_plasmid_qc_701_501_count.csv')
result_dna = get_tables(dna_df)
colAnnon_dna = result_dna[[1]]
dnaCount = result_dna[[2]]
colAnnon_dna<jupyter_output><empty_output><jupyter_text>IMPORTANT ADD IN TISSUE SPECIFIC INFO BELOW  HEK293Tsave_prefix = 'mpraanalyze_HEK293_'sample_ids = c('HEK293_1', 'HEK293_2', 'HEK293_3') Astrocytes (AST)save_prefix = 'mpraanalyze_AST1_'sample_ids = c('AST_1', 'AST_2') H9-D0save_prefix = 'mpraanalyze_H9D0_'sample_ids = c('H9_D0_1A', 'H9_D0_1B', 'H9_D0_2') H9-D2save_prefix = 'mpraanalyze_H9D2_'sample_ids = c('H9_D2_1', 'H9_D2_2', 'H9_D2_3', 'H9_D2_4') H9-D4 (removed H9_D4_3A because of clustering issues)save_prefix = 'mpraanalyze_H9D4_'sample_ids = c('H4_D4_1', 'H9_D4_2', 'H9_D4_3B') H9-D10 (removed H9_D10_2 because of clustering issues)save_prefix = 'mpraanalyze_H9D10_'sample_ids = c('H9_D10_1', 'H9_D10_4', 'H9_D10_3') SLsave_prefix = 'mpraanalyze_SL_'sample_ids = c('SL_1', 'SL_2', 'SL_3A', 'SL_3B', 'SL_4') SLCsave_prefix = 'mpraanalyze_SLC_'sample_ids = c('SLC_1', 'SLC_2', 'SLC_3',  'SLC_4)<jupyter_code>#### SL
save_prefix = 'D_mpraanalyze_barcode_allelic/mpraanalyze_SL_'
sample_ids = c('SL_1', 'SL_2', 'SL_3A', 'SL_3B', 'SL_4')
tissue = 'SL'
rna_df= data.frame()
for (sample_id in sample_ids){
    # read in each lane
    rna_df_rep_L1 = read.csv(paste0('dedup_counts/',sample_id,  '_L1_libpsych_virus_dupumi_count.csv'))
    rna_df_rep_L1 = dplyr::select(rna_df_rep_L1, id, count)
    rna_df_rep_L2 = read.csv(paste0('dedup_counts/',sample_id,  '_L2_libpsych_virus_dupumi_count.csv'))
    rna_df_rep_L2 = dplyr::select(rna_df_rep_L2, id, count)
    
    # merge lanes by summing and add a variable indicating batch
    rna_df_rep = full_join(rna_df_rep_L1,rna_df_rep_L2,by='id')
    rna_df_rep$count.x[is.na(rna_df_rep$count.x)] = 0 
    rna_df_rep$count.y[is.na(rna_df_rep$count.y)] = 0
    rna_df_rep$count= rna_df_rep$count.x +rna_df_rep$count.y
    rna_df_rep$batch= tail(strsplit(sample_id,'_')[[1]],n=1)
    
    #add to overall dataframe
    rna_df = rbind(rna_df, rna_df_rep)
    
}


# get overall dataframe and annotations

result_rna = get_tables_withbatch(rna_df)
colAnnon_rna = result_rna[[1]]
rnaCount = result_rna[[2]]


result_rna_nobatch = get_tables(rna_df)
colAnnon_rna_nobatch = result_rna_nobatch[[1]]
rnaCount_nobatch = result_rna_nobatch[[2]]
snp="chr3_50374568"
rnaCount_nobatch[snp,]
log(sum(as.matrix(rnaCount_nobatch[snp,])[,11:20])/sum(as.matrix(rnaCount_nobatch[snp,])[,1:10]))# log ref/alt
#1.201612 reported 
# filter dna for those found in rna
dnaCount_new = dnaCount %>% 
    rownames_to_column('snpname') %>%
    filter(row.names(dnaCount) %in% row.names(rnaCount))%>%
    column_to_rownames('snpname')

# filter rna for those found in dna
rnaCount_new = rnaCount %>% 
    rownames_to_column('snpname') %>%
    filter(row.names(rnaCount) %in% row.names(dnaCount_new))%>%
    column_to_rownames('snpname')


# filter rna for those found in dna
rnaCount_nobatch_new = rnaCount_nobatch %>% 
    rownames_to_column('snpname') %>%
    filter(row.names(rnaCount_nobatch) %in% row.names(dnaCount_new))%>%
    column_to_rownames('snpname')
rnaCount_nobatch_new = rnaCount_nobatch_new[rownames(dnaCount_new),colnames(dnaCount_new)]
dim(dnaCount)
dim(rnaCount)
dim(dnaCount_new)
dim(rnaCount_new)
dim(rnaCount_nobatch_new)
colAnnon_dna<jupyter_output><empty_output><jupyter_text>MPRAnalyze<jupyter_code># create MpraObject
obj <- MpraObject(dnaCounts = as.matrix(dnaCount_new), rnaCounts = as.matrix(rnaCount_new), 
                  dnaAnnot = colAnnon_dna, rnaAnnot = colAnnon_rna)
# ## Data Normalization
# ## If the library factors are different for the DNA and RNA data, separate 
# ## estimation of these factors is needed. We can also change the estimation 
# ## method (Upper quartile by default)
obj <- estimateDepthFactors(obj, #lib.factor = c("allele", "bc_id"),
                            which.lib = "dna", 
                            depth.estimator = "uq")

obj <- estimateDepthFactors(obj, lib.factor = c("batch"),#lib.factor = c("allele"),
                            which.lib = "rna", 
                            depth.estimator = "uq")
# obj <- estimateDepthFactors(obj, lib.factor = c("allele"),
#                             which.lib = "rna", 
#                             depth.estimator = "uq")<jupyter_output><empty_output><jupyter_text>Quantitative AnalysisThese estimates can then be used to identify and classify SNPs that induce a higher transcription rate (ASE aka Allele specific expression)<jupyter_code># obj <- analyzeQuantification(obj = obj, 
#                               dnaDesign = ~ bc_id + allele,
#                               rnaDesign = ~ allele+batch)<jupyter_output><empty_output><jupyter_text>We can now extract the transcription rate estimates from the model, denoted ‘alpha values’ in the MPRAnalyze model, and use the testing functionality to test for activtiy. extracting alpha values is done with the getAlpha function, that will provide separate values per-factor if a factor is provided. In this case we want a separate alpha estimate by condition:<jupyter_code># ##extract alpha values from the fitted model
# alpha <- getAlpha(obj, by.factor = "allele")
# ##visualize the estimates
# par(bg='white')
# boxplot(alpha)<jupyter_output><empty_output><jupyter_text>Comparative AnalysisMPRAnalyze also supports comparative analyses, in this case: identifying enhancers that are differentially active between conditions. While we can do this indirectly by taking the quantification results and identify enhancers that are active in one condition but not the other, a direct compartive analysis is more sensitive, and allows identification of enhancers that are more or less active, avoiding the binarization of activity. MPRAnalyze also leverages negative controls to estbalish the null differential behavior, thereby correcting any systemic bias that may be present in the data. In terms of syntax, this analysis is done very similarly to quantification, with an additional reduced model that describes the null hypothesis. In this case, the null hypothesis is no differential activtiy between conditions, so the reduced model is an empty model (intercept only) 2212 steps for model takes around ~5-10 min to run<jupyter_code>ptm <- proc.time()
obj <- suppressMessages(analyzeComparative(obj = obj, 
                          dnaDesign = ~ barcode_allelic,
                          rnaDesign = ~ allele, 
                          reducedDesign = ~ 1))
proc.time() - ptm
##
# with the fitted model, we can now test for differential activity, by calling testLrt
res <- testLrt(obj)
head(res)
summary(res)
# for filtering out the snps with < 5 bc
bc_filt_num = 3
rna_df_up_filt = rna_df%>%
    separate(id, c("chr", "loc", "allele_info", "letter", "bc_id"), sep='[.]')%>%
    group_by(chr,loc, allele_info, letter)%>%
    summarise(num_bc = n())%>%
    ungroup()%>%
    filter(num_bc>bc_filt_num) %>%
    mutate(snp_name = str_c(chr, loc, sep = '_'))
snp_names_up_filt = rna_df_up_filt$snp_name
length(snp_names_up_filt)
summary(rna_df_up_filt$num_bc)
ggplot(rna_df_up_filt,aes(x=num_bc))+geom_bar()
# # METHOD 1 dna normalized fold change
# res$normalt = log2(rowSums(rnaCount_nobatch_new[,1:10]+1)/rowSums(dnaCount_new[,1:10]+1))
# res$normref = log2(rowSums(rnaCount_nobatch_new[,11:20]+1)/rowSums(dnaCount_new[,11:20]+1))
# res$logFC = log2(res$normalt/res$normref)

#METHOD 2: just log2 ref/alt in rna
# res$logFC = log2(rowSums(rnaCount_nobatch_new[,1:10]+1)/rowSums(rnaCount_nobatch_new[,11:20]+1))

# # trimming
# res$pval =  pmax(res$pval, 1e-6)
# res$logFC =  pmin(res$logFC, 10)
# res$logFC =  pmax(res$logFC, 1e-6)
dim(res)
res = res%>%
    rownames_to_column()%>%
    filter(rowname %in% snp_names_up_filt)%>%
    mutate(neglog10_p = -log10(pval))%>%
    mutate(neglog10_p_x_logFC = neglog10_p*logFC)%>%
    arrange(desc(neglog10_p_x_logFC))%>%
    drop_na()

# When the hypothesis teseting is simple (two-condition comparison), a fold-change estimate is also available:

## plot log Fold-Change
par(bg='white')
plot(density(res$logFC))
## plot volcano

par(bg='white')
plot(res$logFC, -log10(res$pval))

pdf(paste0(save_prefix, 'volcano.pdf'))
par(bg='white')
plot(res$logFC, -log10(res$pval))

dev.off()
## plot qqplot

# pdf(paste0(save_prefix, 'volcano.pdf'))
par(bg='white')
qq_df = data.frame(emp_pval = sort(res$pval))

qq_df$theoretical_quantile = cume_dist(qq_df$emp_pval)

ggplot(qq_df, aes(x = -log10(theoretical_quantile), y = -log10(emp_pval)))+
    geom_point(size=1) +
    geom_abline(slope=1, intercept = 0, color = scales::brewer_pal(palette = "Set1")(1)[1]) +
  labs(x = "-log10 theoretical p-value", y = "-log10 observed p-value",
       title = paste0(tissue, " QQplot")) +  ylim(0, 30)+
  cowplot::theme_cowplot()
                                                        
ggsave(paste0(save_prefix, 'qqplot.pdf'))
# plot(1:length(-log10(res$pval)), -log10(res$pval))

# dev.off()<jupyter_output>Warning message:
“Removed 1 rows containing missing values (geom_point).”
Saving 7 x 7 in image

Warning message:
“Removed 1 rows containing missing values (geom_point).”<jupyter_text>annotate<jupyter_code>lib_csv = '../../psych_lib_info.csv'
lib_df = read.csv(lib_csv)
lib_df = mutate(lib_df, name_abbr = str_c(Chr, Position, sep='_'))
head(lib_df)
# thres = 0.05/dim(res)[1]
res_filt = res %>%
#     filter(pval < 0.05)
    filter(fdr<0.05)
dim(res)
dim(res_filt)
res_filt = res_filt %>%
    left_join(lib_df, by=c('rowname'='name_abbr'))
res_filt
## plot volcano
# pdf(paste0(save_prefix, 'volcano.pdf'))
# par(bg='white')

ggplot(res, aes(x=logFC, y = -log10(pval))) + geom_point()+
 geom_point(data=res_filt, color='red')+
  labs(x = "log fold change", y = "-log10 p-value",
       title = paste0(tissue, " volcano plot")) + xlim(-5, 5)+ ylim(0,30)+
  cowplot::theme_cowplot()
# plot(res$logFC, -log10(res$pval),title(paste0(tissue, 'Volcano Plot')))
ggsave(paste0(save_prefix, 'volcano.pdf'))

# dev.off()
write.csv(res,paste0(save_prefix, 'comp.csv'))
write.csv(res_filt,paste0(save_prefix, 'comp_filt.csv'))

sort(unique(res_filt$eGene))
sort(unique(res_filt$GeneName))<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup