<p>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.</p>

<p>I assume you would like to get your work published, so going by the papers with IF > 5 such as <a href="http://genesdev.cshlp.org/content/genesdev/32/1/26/F6.expansion.html" rel="nofollow noreferrer">this</a>, <a href="http://genesdev.cshlp.org/content/25/7/742.full" rel="nofollow noreferrer">and this</a>, <a href="https://www.sciencedirect.com/science/article/pii/S1097276515000957" rel="nofollow noreferrer">and this</a>, 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. </p>

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

<pre><code>library(dplyr)
library(tidyr)
library(ggplot2)
library(ggpubr)

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

<p>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. </p>

<pre><code># A tibble: 12 x 3
# Groups:   Transcript [6]
   Transcript   condition value
   &lt;chr&gt;        &lt;chr&gt;     &lt;dbl&gt;
 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
</code></pre>

<p>Then we plot:</p>

<pre><code>p &lt;- 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")))
</code></pre>

<p><a href="https://i.sstatic.net/fZ6MB.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/fZ6MB.png" alt="enter image description here"></a></p>

CC BY-SA 4.0
2020-03-11T15:04:58.263Z
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
kcm
CC BY-SA 4.0
2020-03-11T15:27:19.093Z
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
StupidWolf
CC BY-SA 4.0
2020-03-11T15:27:49.863Z
it's better to spend time on something more interesting like what's in the promoter or SUPER enhancer, than thinking about pausing index
StupidWolf
CC BY-SA 4.0
2020-03-11T16:02:57.047Z
"Warning message:
Computation failed in `stat_signif()`:
not enough 'y' observations " i did get this as im using full data frame
kcm
CC BY-SA 4.0
2020-03-11T16:08:19.563Z
damn i just copy pasted your code and fail to see that...............
kcm
CC BY-SA 4.0
2020-03-11T16:24:00.847Z
so it basically calculates promoter/gene-body the pausing index I do see a significant difference but i read in the paper where they show like the difference in promoter vs gene-body. So can i comment from the above output if the pol-II is high at Promoter or gene-body PTPN6_KO compared to Control condition?
kcm
CC BY-SA 4.0
2020-03-11T16:26:00.100Z
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
StupidWolf
CC BY-SA 4.0
2020-03-11T14:40:03.070Z
2020-03-11T16:06:35.393Z
2020-03-11T16:06:35.393Z
StupidWolf
<p>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.</p>

<p>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? </p>

<p>I'm giving a sample of the output.</p>

<p>Any help or suggestion would be really appreciated. </p>

<pre><code>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))
</code></pre>

<homer>
Pausing index calculation
# Pausing index calculation

(Asked by: username_0 on 2020-03-11)

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))

**Tags:** homer

---

## 1 answers

### Answer by username_1 on 2020-03-11. Score: 1 | Accepted answer

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
                
 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")))

> **Comments:**
> **username_0**: 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  
> **username_1**: 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  
> **username_1**: it's better to spend time on something more interesting like what's in the promoter or SUPER enhancer, than thinking about pausing index  
> **username_0**: "Warning message: Computation failed in `stat_signif()`: not enough 'y' observations " i did get this as im using full data frame  
> **username_0**: damn i just copy pasted your code and fail to see that...............  
> **username_0**: so it basically calculates promoter/gene-body the pausing index I do see a significant difference but i read in the paper where they show like the difference in promoter vs gene-body. So can i comment from the above output if the pol-II is high at Promoter or gene-body PTPN6_KO compared to Control condition?  
> **username_1**: 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  

---

