2
$\begingroup$

I have Pol_II data for Ctrl vs Knockout condition with two replicates each. I want to calculate the pausing index for the Ctrl vs KO condition.

Here I have given tag directory as input to Homer.Now what I see is as such,for each biological replicate it calculates the pausing ratio.How can I calculate as a whole like 2 Ctrl vs 2 KO replicate?

I'm giving a sample of the output.

Any help or suggestion would be really appreciated.

dput(head(CTRL_PTPEM_Pol_II_pausing))
structure(list(`Transcript/RepeatID (cmd=analyzeRepeats.pl rna hg38 -strand both -count pausing -condenseGenes -d Control_Pol_II_new_set/ Control_Pol_II_set2/ PTPN6_KO_Pol_II_PTPEN/ PTPN6_KO_Pol_II_set2/ -rpkm)` = c("NM_173803", 
"NM_001301039", "NR_132775", "NM_022114", "NM_001011666", "NR_160728"
), chr = c("chr16", "chr1", "chr3", "chr1", "chr7", "chr19"), 
    start = c(15395704, 150272234, 13618331, 3069153, 28686083, 
    22527485), end = c(15400754, 150277284, 13623381, 3074203, 
    28691133, 22532535), strand = c("+", "+", "+", "+", "+", 
    "-"), Length = c(250, 250, 250, 250, 250, 250), Copies = c(2, 
    1, 1, 1, 1, 1), `Annotation/Divergence` = c("MPV17L|M-LPH|MLPH1|MLPH2|MPV17L1|-|16p13.11|protein-coding", 
    "C1orf54|-|-|1q21.2|protein-coding", "SNORA93|-|-|3p25.1|snoRNA", 
    "PRDM16|CMD1LL|KMT8F|LVNC8|MEL1|PFM13|-|1p36.32|protein-coding", 
    "CREB5|CRE-BPA|CREB-5|CREBPA|-|7p15.1|protein-coding", "LOC105376917|-|-|19p12|ncRNA"
    ), `Control_Pol_II_new_set/ Pausing Ratio` = c(0.982, 1.832, 
    1.114, 1.247, 3.476, 0.632), `Control_Pol_II_new_set/ Promoter Reads (-50,200)` = c(0.125, 
    0.748, 0.249, 0.873, 0.499, 0), `Control_Pol_II_new_set/ GeneBody Reads (200,5000)` = c(0.13, 
    0.338, 0.208, 0.669, 0.032, 0.091), `Control_Pol_II_set2/ Pausing Ratio` = c(0.686, 
    1.131, 1.054, 0.816, 2.205, 1.005), `Control_Pol_II_set2/ Promoter Reads (-50,200)` = c(0.1, 
    0.7, 0.1, 0.3, 0.3, 0.1), `Control_Pol_II_set2/ GeneBody Reads (200,5000)` = c(0.203, 
    0.604, 0.089, 0.396, 0.068, 0.099), `PTPN6_KO_Pol_II_PTPEN/ Pausing Ratio` = c(1.08, 
    1.362, 1.108, 1.224, 1.943, 1.2), `PTPN6_KO_Pol_II_PTPEN/ Promoter Reads (-50,200)` = c(0.151, 
    0.602, 0.151, 0.301, 0.452, 0.151), `PTPN6_KO_Pol_II_PTPEN/ GeneBody Reads (200,5000)` = c(0.125, 
    0.392, 0.118, 0.212, 0.141, 0.094), `PTPN6_KO_Pol_II_set2/ Pausing Ratio` = c(1.652, 
    0.909, 1.551, 1.395, 1.774, 0.584), `PTPN6_KO_Pol_II_set2/ Promoter Reads (-50,200)` = c(0.429, 
    0.536, 0.429, 0.536, 0.321, 0.107), `PTPN6_KO_Pol_II_set2/ GeneBody Reads (200,5000)` = c(0.206, 
    0.603, 0.229, 0.346, 0.123, 0.279)), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -6L))
$\endgroup$
1
$\begingroup$

but I would say the pausing index is something that is invented, commonly used, but with no clear indication of what is the expected distribution or how it can be subjected to test like you've mentioned.

I assume you would like to get your work published, so going by the papers with IF > 5 such as this, and this, and this, the clear strategy is use a boxplot to show the difference between conditions, and apply a statistical test such as mann whitney u or t.test to say, there is a difference.

Let's try this, sorry not very good in R,using your dataframe which i named x:

library(dplyr)
library(tidyr)
library(ggplot2)
library(ggpubr)

df = x %>% select(contains("pausing") | contains("Transcript")) %>% 
rename_at(vars(1),substr,1,10) %>% 
pivot_longer(-Transcript,names_to = "condition") %>% 
mutate(condition=case_when(
grepl("Control",condition) ~ "Control",
grepl("PTPN6",condition) ~ "PTPN6_KO"
)) %>%
group_by(Transcript,condition) %>% summarize_all(mean) 

In the data frame above, I have averaged the pausing index between replicates (not very good but then again, most papers have only 1 replicate), and group them by condition.

# A tibble: 12 x 3
# Groups:   Transcript [6]
   Transcript   condition value
   <chr>        <chr>     <dbl>
 1 NM_001011666 Control   2.84 
 2 NM_001011666 PTPN6_KO  1.86 
 3 NM_001301039 Control   1.48 
 4 NM_001301039 PTPN6_KO  1.14 
 5 NM_022114    Control   1.03 
 6 NM_022114    PTPN6_KO  1.31 
 7 NM_173803    Control   0.834
 8 NM_173803    PTPN6_KO  1.37 
 9 NR_132775    Control   1.08 
10 NR_132775    PTPN6_KO  1.33 
11 NR_160728    Control   0.818
12 NR_160728    PTPN6_KO  0.892

Then we plot:

p <- ggboxplot(df, x = "condition", y = "value",
                color = "condition", palette =c("#c72e29","#016392","#be9c2e","#098154"),
                add = "jitter")
p+stat_compare_means(comparisons=list( c("Control", "PTPN6_KO")))

enter image description here

$\endgroup$
7
  • $\begingroup$ your response are life saver let me go through the links and will update you. One doubt can you just take average of the replicates thought of doing that but was not sure if that is right thing because homer might have some model to calculate it $\endgroup$
    – kcm
    Mar 11 '20 at 15:04
  • 1
    $\begingroup$ you're welcome. well let me know if you find a way. if you just want to make a statement like it pauses more in condition a vs b, i think the above is definitely good enough $\endgroup$
    – StupidWolf
    Mar 11 '20 at 15:27
  • 1
    $\begingroup$ it's better to spend time on something more interesting like what's in the promoter or SUPER enhancer, than thinking about pausing index $\endgroup$
    – StupidWolf
    Mar 11 '20 at 15:27
  • $\begingroup$ "Warning message: Computation failed in stat_signif(): not enough 'y' observations " i did get this as im using full data frame $\endgroup$
    – kcm
    Mar 11 '20 at 16:02
  • 1
    $\begingroup$ you can make a meta-gene plot, combined with this pausing index above, will tell you, most likely you have more pol-II at stuck at the promoter... i think it's pretty obvious, your phenotype $\endgroup$
    – StupidWolf
    Mar 11 '20 at 16:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.