text:
GC normalization ATAC seq data
Is it really required to do GC normalization ATAC seq data ?
One of the paper where they have ATAC seq data did the normaisation of GC bias after peak calling. When i looked for the library and the reason why GC bias i found this
Could it be true that genes with higher GC content are higher expressed? It has been suggested that genes that are either extremely high or extremely low expressed are under some form of selection leading to “extreme” GC content. What CQN does, is making the effect of GC content comparable across samples, and we show in [1] that this leads to improved inference. It also flattens the effect of GC content on gene expression, but we believe this is better than having the effect of GC content depend on the sample.
https://www.bioconductor.org/packages/release/bioc/vignettes/cqn/inst/doc/cqn.pdf

Now that scenario is for RNA seq would it be logical to do the same for ATAC seq?

" I personally do not do it. It is not established from what I know and I do not see the advantage. " i have the doubt but not the knowledge  if i should do or not ! ..i came across paper where they did that but not sure why ..but now im clear how to proceed thank you and @Devon Ryan for the answer

The only reason to normalize for GC content in RNA-seq is if it differs notably between samples/groups. If that's not the case and you aren't trying to compare genes withing samples then you have no reason to try to account for GC content. The same goes for ATAC-seq, though there the danger with trying to normalize for GC content is that you then mask changes in global accessibility.
In short, normalize for things that should cause problems with your analysis and nothing else.

Now my confusion is cleared


text:
the model matrix is not full rank : This is a classic question which a biologist face without clear understanding of the model design
Saw this answer biostar. Tried to make my metadata as such but still the
"Error in checkFullRank(modelMatrix) : "

This is my coldata
dput(coldata)
structure(list(Group = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L,
6L, 6L, 6L, 6L), .Label = c("HSC", "Blast", "CMP", "GMP", "LSC",
"Mono"), class = "factor"), Mutation = structure(c(10L, 4L, 3L,
5L, 2L, 3L, 9L, 11L, 1L, 1L, 7L, 7L, 10L, 7L, 8L, 2L, 10L, 10L,
10L, 10L, 2L, 2L, 7L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 2L, 2L, 5L, 2L,
9L, 1L, 1L, 8L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("DNMT3A-R882H",
"FLT3-ITD", "IDH2-R140Q", "KRAS-G13D", "MAU2-Q98H", "No", "NPM1-L287ins(TCTG)",
"NRAS-G13D", "STAG2-R614", "TET2-E1357", "TET2-R550"), class = "factor"),
    Satus = structure(c(2L, 5L, 4L, 6L, 3L, 4L, 9L, 10L, 1L,
    1L, 7L, 7L, 2L, 7L, 8L, 3L, 2L, 2L, 2L, 2L, 3L, 3L, 7L, 16L,
    16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L,
    16L, 16L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 12L, 12L, 13L,
    12L, 9L, 11L, 11L, 14L, 15L, 15L, 15L, 15L, 15L, 15L), .Label = c("Blast-DNMT3A-R882H",
    "Blast-E1357", "Blast-FLT-ITD", "Blast-IDH2-R140Q", "Blast-KRAS-G13D",
    "Blast-MAU2-Q98H", "Blast-NPM1-L287ins(TCTG)", "Blast-NRAS-G13D",
    "Blast-STAG2-R614", "Blast-TET2-R550", "LSC-FLT-DNMT3A-R882H",
    "LSC-FLT-ITD", "LSC-MAU2-Q98H", "LSC-NRAS-G13D", "Mature-Normal",
    "Progenitor-Normal", "Stem-Normal"), class = "factor")), row.names = c("Blast1",
"Blast10", "Blast11", "Blast12", "Blast13", "Blast14", "Blast15",
"Blast16", "Blast17", "Blast18", "Blast19", "Blast20", "Blast2",
"Blast21", "Blast22", "Blast23", "Blast3", "Blast4", "Blast5",
"Blast6", "Blast7", "Blast8", "Blast9", "CMP1", "CMP2", "CMP3",
"CMP4", "CMP5", "CMP6", "CMP7", "CMP8", "GMP1", "GMP2", "GMP3",
"GMP4", "GMP5", "GMP6", "GMP7", "HSC1", "HSC2", "HSC3", "HSC4",
"HSC5", "HSC6", "HSC7", "LSC1", "LSC2", "LSC3", "LSC4", "LSC5",
"LSC6", "LSC7", "LSC8", "Mono1", "Mono2", "Mono3", "Mono4", "Mono5",
"Mono6"), class = "data.frame")

dds <- DESeqDataSetFromMatrix(countData=df2, colData=coldata, design= ~ Group + Mutation)

This is the design I'm trying to incorporate. What i get is my design is nested in one of the column by looking at numerous example.
What i would like to do is I'm trying to run a LRT test on this dataset which works fine if i give "Group" as my design and then i run this
dds_lrt <- DESeq(dds, test="LRT", reduced = ~ 1)

which gives me diffrenitaly expressed across multiple groups. Now there are two question

Is it possible to include interaction such as in my case Group+Mutation and do the lrt test?
Or the above is conceptually wrong to do ?

what is the hypothesis you are trying to test?

I would like to know how these specific mutation affect gene expression across the cell types with respect to reference which is my "HSC" as its a stem cell . But from @ATpoint answer it seems I need to subset the data more over what i tired above perhaps is a wrong approach if i get it as there are multiple types of mutation. It would be more useful if i do pairwise comparison for specific cell types i guess

@StupidWolf for example I would like to see what set of genes which changes or are significant when i traverse from HSC to down the lineage tree such as CMP,GMP,Monocyte, but I have also included LSC(Leukemic stem cells) and Blast. But not sure if this is right comparison, i did normalize and took the data for clustering that is fine ,but hypothesis testing across heterogeneous group which are diseased with mutation would that be correct?

@StupidWolf this answer let me go through again https://bioinformatics.stackexchange.com/questions/13627/how-to-design-deseq2-lrt-model-with-individuals-nested-in-2-levels

Looking at your design with ATpoint's table, there's a major issue. some of your mutations have n=1 in some of the interactions, for example TET2-R550 or NRAS-G13D. or NPM1-L287ins(TCTG) is exclusive for blast. So it's hard for you to make these comparisons say HSC.. when it is so incomplete

I think you basically need to think about how to design this

"there's a major issue. some of your mutations have n=1 in some of the interactions, for example TET2-R550 or NRAS-G13D. or NPM1-L287ins(TCTG) is exclusive for blast. So it's hard for you to make these comparisons say HSC" now i get it at least i need a replicate ..may be i can make another categorical variable

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Open chromatin region : nucleosome , mono-nucleosome , di-nucleosome
I m following this tutorial for my ATAC seq analysis. Here in of the section they show profile of open chromatin region plots.
**Plotting insert sizes with Greenleaf open, mono- and di-nucleosome profiles
This looks very similar to the image from the Greenleaf paper.

We can now annotate our nucleosome free (< 100bp), mono-nucleosome (180bp-247bp) and di-nucleosome (315-437) as in the Greenleaf study.**

Now Since i have done the differential accessibility analysis. So from that information how do i find out region which are "nucleosome free (< 100bp), mono-nucleosome (180bp-247bp) and di-nucleosome (315-437) ".
My objective

I would like to see this "nucleosome free (< 100bp), mono-nucleosome (180bp-247bp) and di-nucleosome (315-437) " from my differential accessible region which was done from readcount and then I ran deseq2.
My goal is not to plot the way they have done rather find those regions mentioned above and see it between different conditions.

Now based on the tutorial of that section where they do for each bam file.
Do i need to run the analysis on each bam file or each condition or Can I find out the above from my deseq2 analysis which has regions and those regions I have annotated it .

It's impossible to find nucleosome free and such regions without a BAM file. Your only options are:

Filter the BAM files to contain only the size-range of interest and rerun the counting and statistics.
Find or write a counting program that can track fragment sizes.

BTW, those size ranges could use some adjustments, you're throwing away all fragments between 100 and 180 bases using their ranges.

Thank you devon . "BTW, those size ranges could use some adjustments, you're throwing away all fragments between 100 and 180 bases using their ranges." how to do it add some region?

Looks at a size distribution plot and adjust accordingly.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
ATAC seq : Plotting insert sizes
My tutorial link
library(GenomicAlignments)

atacReads <- readGAlignmentPairs("MACS2_peak_call/ALL_ATAC_data/HSC1.bam", param = ScanBamParam(mapqFilter = 1,
            flag = scanBamFlag(isPaired = TRUE, isProperPair = TRUE), what = c("qname",
        "mapq", "isize"), which = GRanges("chr20", IRanges(1, 63025520))))
# length(atacReads)
atacReads

atacReads_read1 <- GenomicAlignments::first(atacReads)
insertSizes <- abs(elementMetadata(atacReads_read1)$isize)
head(insertSizes)

library(magrittr)
library(dplyr)
library(ggplot2)
fragLenPlot <- table(insertSizes) %>% data.frame %>% dplyr::rename(InsertSize = insertSizes,
    Count = Freq) %>% mutate(InsertSize = as.numeric(as.vector(InsertSize)),
 Count = as.numeric(as.vector(Count))) %>% ggplot(aes(x = InsertSize, y = Count)) +
  geom_line()

fragLenPlot + theme_bw()

fragLenPlot + scale_y_continuous(trans = "log2") + theme_bw()

How to i specify the Grange to span the total chromosome. Now I use only "chr20" as shown in the example.

sorry don't quite get the question. is it about knowing the chromosome length from the bam file, or it's something to do with the plotting

Now Im too bit confused ,after reading paper where they don't mention which chromosome they taking and plotting.  It would be essentially this one " is it about knowing the chromosome length from the bam file" can i like  "chr1-chr20" and see the insert size distribution . I was not able to figure out how to define more that one chromosome in that Grange argument //

@StupidWolf for example in this paper https://www.nature.com/articles/s41597-019-0071-0#Sec1  Fig2(b).

You can exclude the which statement and get the entirety of all of the chromosomes. Note however that R is a very bad platform to use for this, since you end up reading everything into memory (I hope the BAM file isn't huge).
Note that you're seriously reinventing the wheel. There are a number of packages that will directly produce the plot you want, such as bamPEFragmentSize from deepTools.

"bamPEFragmentSize" will see this

" you're seriously reinventing the wheel." yes since I didn;t read about your tool "bamPEFragmentSize" I used it , Works absolutely fine

is there any way by using this "bamPEFragmentSize" i can get smooth curve ?

Save the values to a file and replot however you like.

okay i will see the option to save the file ....right now i ran this "bamPEFragmentSize -hist fragmentSize.png --maxFragmentLength 100 -b HSC1.bam Mono1.bam LSC1.bam Blast1.bam --numberOfProcessors 10 --samplesLabel HSC Mono LSC Blast
"


text:
Nanopore data clarification regarding merging of samples
We had same set of sample ran on nanopore platform twice.
The issue on hand is.
First run was done it was nearly done but the power went off.
Second run was done completely with enough sequencing depth.
Now the question I ran the analysis with the second set which was complete that set of result I have. But so far I have't ran the first part.
Question1 : if I run the first part do i get major differences? in the result
Question2 : would it be advisable to merge both run since they are now can be called as twi different batches plus would be bring drastic changes in the output ? And in terms of tool for merging I googled which shows poretools.

It is entirely unclear which type of analysis you aim to do. You probably should be more specific when asking questions.
That said, given that the run was going perfectly fine until a power issue stopped it, all data produced should be fine and you can indeed combine those runs.
You probably shouldn't bother with poretools. If you are using fastq files to start your analysis you can just cat the basecalled files together. Also after alignment, you can merge the bam files. If your analysis requires the fast5 files then you can probably just put these in the same folder, but you should tell us which analysis you are doing.

So we have gut mircobiome sample ..between two groups. For classification Im using kraken2 and the output is used for downstream analysis using phyloseq


text:
ATAC seq density calculation
This paper Fig2b

b Changes of chromatin accessibility in AMD (n = 14) relative to
normal (n = 11) in all retina samples. Each dot represents one
ATAC-Seq peak. Blue line in the left panel indicates average fold
changes of peaks with the same ATAC-Seq intensity. The percentage of
reduced peaks is shown under the density curve in the right panel.

How is the density calculated or what is being shown there? average ATAC seq signal or Fold change?

They are not referring to density of reads. It is the density of the fold change, the y variable in the plot to left of Fig2b. It tells you the distribution of the foldchange and 92% of the fold change is <0.
See this page for more details on density and histograms
Suppose your data is something like this (purely making this up):
set.seed(555)
fc = rnorm(1000,-0.5,2)
averATAC = -0.5*(fc - 0.5)^2 + 0.2*fc + 3  + rnorm(100,25,10)

You fit a density to it and you get the plot:
par(mfrow=c(1,2))
plot(averATAC,fc)
plot(dens$y,dens$x,xlab = "density",ylab="fc")

your answers are part of my phd work I mean those example which helps me to express my solution

how is the percentage being calculated from the density?

I think it's an overkill, but you can do sum(dens$y[dens$x<0])/sum(den$y)

"overkill" i will go for it

Just count the number of reads overlapping an ATAC-seq peak and divide it by the length of the peak gives you the read density. Each dot is a fold change, but to calculate fold change, average density from all the samples are taken, that is what the "average" is referring to.

my query is about the 92.3 % density what does it represent ? The fold change or the average ?

I see, it's the proportion of fold change falls below 0 on the y axis of the plot on the left

so is that the negative fold change they plotted?

it is a density plot of fold changes, not the density of the reads.


text:
Kruskal–Wallis and wilcox.test choice of test for expression data
I have bit of confusion regarding how ggpubr has implemented the KW and wilcox.test. By default usage of ggpubr reports p value for wilcox.test() but when it comes to Kruskal–Walis it is different .
But as i have to compare more that two groups im opting for Kruskal–Wallis so now the confusion is i get a p value which is significant but if i do the test independently then the p value is not significant not sure why is that
I have asked similar question old question
My data
dput(df1)
structure(list(mHSC = c(7.7559725, 3.307901197229, 7.91049997141824,
10.9659075, 5.66225942792448, 9.40594075, 8.03212169895329, 10.929870369027,
9.80897606357472, 4.29298594445745, 8.58579242742772, 12.0936123426199,
8.61371460307564, 11.0174749823542, 10.5147718423478, 5.53565725,
8.12612925, 10.6881344317581, 8.48826962785442, 5.52294757621421,
8.38046925, 9.44552809599148, 9.19142493299942, 7.03072998196466,
8.88679975, 0.3126697975, 10.4907798679713, 7.864884, 12.7182883447252,
6.3180225, 5.87640020441158, 7.87333318550091, 8.15221192095875,
4.51086624773858, 5.957332, 6.90057575, 5.64708909092112, 10.7124580115802,
9.82466018308675, 8.49078, 8.63046025, 6.96721675, 7.57871225,
8.815752, 6.72449933154642, 10.7951750189282, 8.01780331997795,
6.07469990935439, 8.31037653431384, 9.60791091202699, 7.11738304581177,
6.27364283988744, 10.8712275, 4.76963356477331, 8.22113922937348,
4.68111125, 7.81847825, 8.4123615, 7.1236546846872, 5.5773485,
10.72958, 7.716703, 7.606539, 9.3070835272045, 7.75540042667672,
9.8687305, 8.98201221147182, 8.93182675, 10.868215, 10.8990875,
10.9377104128057, 8.60048104950787, 7.81762687302046, 9.24377416902964,
6.57057975, 7.52549392459474, 5.9971085, 10.7335671888057), mCMP = c(8.08511525,
3.15081874988432, 8.21803244045432, 12.7012625, 6.95234868015311,
10.24927775, 8.02201688082928, 11.1415993774267, 8.6576364168622,
4.33979841956685, 9.20120340882164, 12.6360469440726, 8.82394480426142,
11.7614369660881, 10.518813583537, 5.9224525, 8.3939755, 10.9519826793289,
9.07562352307607, 4.68003855905254, 8.7914695, 9.35721245249627,
9.52072790562384, 7.52153767036027, 9.37291475, 5.323077, 11.316187036788,
8.08079325, 12.9250246147509, 6.3281185, 6.36628784788811, 8.26354221670866,
8.80757801591239, 5.08139074215397, 6.6258195, 7.56676875, 6.30618317691088,
12.8268812686057, 10.4361513118013, 9.12211125, 8.94618675, 7.25775925,
8.7313635, 9.8282665, 7.10318403743241, 10.9016853993601, 8.75786307748107,
6.2613305781586, 7.24503928286349, 10.1900260653741, 6.04352467875968,
5.41919427988363, 11.689205, 5.2447156735349, 7.81023802636728,
5.484589, 7.9321635, 8.8598865, 7.26208235921707, 11.7873, 11.256575,
8.574748, 8.3324525, 9.69726624009281, 8.51491553090231, 10.792235,
9.97113940283311, 10.08911925, 11.867595, 11.8185525, 11.9422571106972,
9.19311951548278, 8.61818925695215, 9.34701182556221, 7.063838,
7.41275116679947, 6.84360325, 10.478012600725), mGMP = c(8.07654075,
3.12607874540186, 7.91384610943697, 12.9717825, 6.533351883655,
10.394515, 8.26093661667729, 11.2346407524235, 11.1019372952288,
4.34824706553879, 9.19896153987714, 12.6475029547774, 8.53466968298702,
11.6750087661439, 10.7705094222618, 5.60385575, 8.32733375, 10.6929416341222,
8.85393492258581, 4.40609754458633, 9.2379455, 9.43854664655775,
9.66597253836355, 7.99640783151295, 9.60382225, 0.40676415, 11.3480173068737,
8.2298565, 13.3988341481968, 5.2305385, 6.12320931339542, 8.16212407932226,
8.98321413854527, 5.21894298832466, 6.88443175, 7.9935035, 6.66952804861909,
11.3240966560241, 10.795029283416, 8.830129, 8.58251875, 7.6113605,
8.508246, 9.63661975, 6.87344684043907, 11.2237814311451, 9.43948359781618,
6.10350570815157, 8.37769527299503, 9.98272594007505, 7.15195399920282,
5.48458261850125, 11.8302525, 5.66643624946888, 8.36658108147322,
5.9928725, 8.56603575, 9.670526, 7.70714893990735, 9.552538,
11.5112725, 8.99551975, 9.504391, 9.76913665769114, 8.58112436662299,
11.2680775, 10.5535730280208, 10.226635, 11.9211375, 12.0610625,
12.3249064903557, 9.52625419135532, 8.66295412469679, 10.0469451814325,
7.700204, 7.44384030300093, 7.40615025, 10.222430040369), mMono = c(6.86074675,
0.899186661852538, 5.72102113412742, 11.0054, 5.22275722758541,
8.56569375, 7.25543228485158, 9.37557140457706, 7.92766781596895,
1.22333714157354, 7.54589243251616, 11.0881258432345, 7.91303390155615,
10.3482735440588, 7.2731491579751, 3.10091425, 6.91085225, 9.03107467062931,
7.91879733419338, 2.71978845710368, 8.31396425, 6.89137159431571,
8.68813944060236, 6.65073272892769, 8.41960525, 0.283334725,
10.0529313221945, 6.7606855, 11.6800576702458, 3.6173355, 3.81146596960356,
6.93219836413735, 7.92727013246883, 4.09314670396697, 5.78634575,
6.55346475, 5.05937426539404, 10.0514401266737, 9.21573912287031,
7.80863325, 7.80450125, 5.96109, 6.07350675, 7.14264825, 5.20991063841094,
8.9541923959891, 5.48494732937368, 3.46755807112085, 5.91011468781522,
6.19799161521768, 4.09986325310833, 1.58331376348435, 11.021745,
4.62557956951436, 7.48457597125806, 4.9819695, 7.567706, 8.547291,
7.13353425765299, 9.4336955, 10.6991475, 7.3555725, 8.32300975,
9.14731426222114, 8.2750379299376, 10.85821, 9.63317457338075,
9.6021585, 11.4682025, 11.4990425, 11.4153866521691, 9.02074441245752,
7.96069342402473, 9.45598020603536, 6.462563, 6.78696488936524,
5.50898575, 8.79251593030914)), class = "data.frame", row.names = c(NA,
-78L))

Code for the first figure with wilcox test and KW as argument
  head(df1)

ex <- melt(df1, id.vars=c(NULL))

head(ex)
my_comparisons <- list( c("mHSC","mCMP"),c("mCMP","mGMP"),c("mGMP","mMono"),c("mHSC","mMono"),c("mHSC","mGMP"))

head(ex)

ggboxplot(ex, x = "variable", y = "value",ggtheme = theme_bw(base_size = 30),
          color = "variable", palette = "jco", add = "jitter")+
  stat_compare_means(comparisons = my_comparisons,method = "kruskal.test")+ # Add pairwise comparisons p-value
  stat_compare_means(label.y = 50)     # Add global p-value

New figure
The one im not able to get is in the first figure even with Wilcox and KW as argument the global p value for KW doesn't change where as the KW p value for comparison gives p value which is not significant .
Am i doing it correctly ?

The lower graph doesn't look good. The "Wilcoxan values" are all identical and I can't see how if all pairwise comparisons were 0.46 that would result in a KW of <<0.001. The way the values hold together in the upper graph looks fine.
Just to summarise

KW = 1 way ANOVA ... non-parametric version
Wilcoxan = T test ... non-parametric version

Looking at your code you specify the pairwise comparisons but then in the ggplot2 code you state method = "kruskal.test" ... are you sure that isn't method = "wilcoxan.test"?
KW will want to perform a global comparison across multiple single columns of data and this might have resulted in some sort of weird mean value being generated, because it will just want to output a single value (my bet is that is why all values are identical). Wilcoxan will be expecting to a request to perform pairwise tests.

Okay my advice is to perform an independent Wilcoxan test outside of ggpubr
wtest <- wilcox.test(mCMP ~ mHSC, data = df1, paired = TRUE)
wtest

Output
    Wilcoxon signed rank test
data:  weight by group
.....
wtest$p.value

Output
p value
I would then simply remove the calculation from with the boxplot and present the p values independently. The only other possibility is whether the data is paired or not and you may have changed that setting between calculations.

"method = "kruskal.test"  are you sure that isn't method = "wilcoxan.test" i read in the ggpubr you can specify the method with argument "kruskal.test" ,by default it does wilcox test."Wouldn't you simply perform pairwise T-test" well the data doesn;t come from normal distribution so i didn't perform t-test . KW test because i have more than two comparison ."KW will want to perform a global comparison across multiple single columns of data,this might have resulted in some sort of weird mean value being generated, we will get reviewers comment why wilcox test because of the number of sample.

Trouble shooting described below. This isn't a biggy.

will use it as you mentioned but a bit of doubt "imply remove the calculation from with the boxplot and present the p values independently" so this is acceptable in publications i have seen though. But if i get it its better to do a pairwise Wilcox rather KW global test even though im going for multiple comparison

right i ran into this error "Error in wilcox.test.formula((mHSC ~ mCMP), data = df1, paired = TRUE) :
grouping factor must have exactly 2 levels"

this works "wtest <- wilcox.test(df1$mHSC , df1$mCMP, paired = TRUE)"

Cool, and now you've mentioned it, this makes sense using dataframe.column format.


text:
Regular expression struggle
These are my files
SO_8543_IP_3_ablt_R1.fq.gz
SO_8543_IP_4_ablt_R1.fq.gz
SO_8543_INPUT_1_ablt_R1.fq.gz  SO_8543_IP_5_ablt_R1.fq.gz
SO_8543_INPUT_2_ablt_R1.fq.gz  SO_8543_IP_6_ablt_R1.fq.gz
SO_8543_INPUT_3_ablt_R1.fq.gz  SO_8543_IP_7_ablt_R1.fq.gz
SO_8543_INPUT_4_ablt_R1.fq.gz  SO_8543_IP_8_ablt_R1.fq.gz
SO_8543_IP_1_ablt_R1.fq.gz
SO_8543_IP_2_ablt_R1.fq.gz

Im running the below line
for i in $(ls *.fq*.gz | cut -f1-5 -d"_" | uniq);do echo $i; done

bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38  -1 ${i}_R1.fq.gz | samtools view -bS - >${i%}.bam;

done

So if i use run this
for i in $(ls *.fq*.gz | cut -f1-5 -d"_" | uniq);do echo $i; done

I get files as such
SO_8543_INPUT_1_ablt
SO_8543_INPUT_2_ablt
SO_8543_INPUT_3_ablt
SO_8543_INPUT_4_ablt
SO_8543_IP_1_ablt
SO_8543_IP_2_ablt
SO_8543_IP_3_ablt
SO_8543_IP_4_ablt
SO_8543_IP_5_ablt
SO_8543_IP_6_ablt
SO_8543_IP_7_ablt
SO_8543_IP_8_ablt

The part of the file name that is consistent is _R1.fq.gz which goes after the input argument. As these are single end files so only 1 im giving one end as input.
Error message
***
Error: Must specify at least one read input with -U/-1/-2
(ERR): bowtie2-align exited with value 1
./Chip_SEQ.sh: line 25: syntax error near unexpected token `done'
./Chip_SEQ.sh: line 25: `done'

it looks like its not able to recognize as input .Any help or suggestion would be highly appreciated

You have an extra done, instead you want:
for i in $(ls *.fq*.gz | cut -f1-5 -d"_" | uniq); do
    echo $i
    bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 -1 ${i}_R1.fq.gz | samtools view -bS - >${i}.bam
done

As an aside, I'm not a big fan of piping ls to cut and then uniq. Why not instead:
for f in *.fq.gz; do
    sampleName=${f%%_R1.fq.gz}
    echo $sampleName
    bowtie2 -q ...
done

That will do the same thing, but (A) it's clearer what's going on and (B) it can be reused if your sample names change. Also, stop writing bash scripts and use snakemake or something like that instead, it scales much better and is less error-prone.

now after looking at my script number of times i figured out the "done" is extra .let me try your code. My shell scripting is very rudimentary so i use which i know is safe to save my self from screwing up .

i ran your first code the but i still get this error "Error: Must specify at least one read input with -U/-1/-2 (ERR): bowtie2-align exited with value 1" its  not able to recognize the input file, not sure what's wrong .

Try -U instead of -1 in the bowtie2 command.

yes that was the mistake i was doing thank you very much

@krushnachChandra please read this link to understand why you should basically never try to parse the output of ls.

@terdon I'd not say "never parse ls", I'd say "Don't parse ls unless you know exactly what you're doing, especially the filenames involved". I do a lot of for f in $(/bin/ls -1 <glob> | downstream_operations) because I'm 100% certain that each value of $f plays well with the other parts in the code.

@RamRS well, parsing ls is asking for trouble and just unnecessary. for f in <glob> does exactly the same thing but without calling an external process and without running any risk of failure because of file names. So why use a riskier and less efficient approach when it doesn't provide any benefit at all?

It's really useful when the $() feeding the loop involves a uniq, for example. It can't always be avoided. Whenever possible, use a glob. But don't have a hard constraint against $(/bin/ls | downstream_stuff)

thats a life saver link


text:
Intersecting multiple data frames
I have to intersect multiple data frame to find common elements. For visualization purpose i can use upset Plot but it doesn;t return a object containing the intersection object .I tried jvenn but at max it can take six sets to intersect .
Im able to do the interesction but i would like to know if i can do it more efficiently ,which apparently not sure how to do this is my code as of now where i read each file as data frame, covert them to list and pass them to do intersection .
NewList <- split(GO0000122, f = seq(nrow(GO0000122)))
NewList1 <- split(GO0006351,f=seq(nrow(GO0006351)))
NewList2 <- split(GO0006355,f=seq(nrow(GO0006355)))
NewList3 <- split(GO0006357,f=seq(nrow(GO0006357)))
NewList4 <- split(GO0006366,f=seq(nrow(GO0006366)))
NewList5 <- split(GO0030154,f=seq(nrow(GO0030154)))
NewList6 <- split(GO0045892,f=seq(nrow(GO0045892)))
NewList8 <- split(GO0045893,f=seq(nrow(GO0045893)))
NewList7 <- split(GO0045944,f=seq(nrow(GO0045944)))

a <- intersect_all(NewList,NewList1,NewList2,NewList3,NewList4,NewList5,NewList6,NewList7,NewList8)

b <- do.call(rbind.data.frame, a)

Can i read the files and pass them into list and do the same that would help me from doing lot of things manually
I can read the files from the directory where i have my files but not sure how to pass them to make it as list.
filenames <- gsub("\\.csv$","", list.files(pattern="\\.csv$"))
for(i in filenames){
  assign(i, read.csv(paste(i, ".csv", sep="")))
}

The object filenames is character here .It just returns me list of files.
So how to make my question bit more clear , i would like to read the list of file from the directory and then pass these to perform intersection. Any suggestion or help would be really appreciated .

I've been cleaning some of your post for punctuation, but this one I am not sure how to change it. Could you please [edit] and correct the comas and i to I and other punctuation errors? Also, if you could clarify what do you mean by intersection here I think that your problem might be easier to solve if you explain what is your goal with this (there might be an out of the box solution)

oh sorry typing partly from phone adds to the problem. My goal is to get set of common elements from multiple data frame where their is only single column consist of genes.So as of now i read each file,then pass those as list, then i do intersect and then turn it into dataframe. So now each time i have to read files if i have 12 files i have to read 12 times .Is there a way to put it in a loop and do the same what i have done above.

I think it would be it easier to work with character vectors than lists (it will at least make it easier to see them). But what have to do the files with the intersection, if any first you read and then you intersect, not the other way around or I am missing something here?

it would be easier if i add my files which im using  https://drive.google.com/open?id=1laieKec1sr4rWa4SIhUvHcLSerrWf4OC

Easier for you but not for me. I don't understand what have to do this files with your problem

okay may be its too confusing but can you let me know if i can how can i read all the csv files in a loop and make them into data-frame that would be easier for me to do the rest

This question is out of the scope of this site. You can find your solution at StackOverflow, specifically this question

sure this helps

Change dataframe to load all the files into a list:
wd <- "C:/all_your_dataframes_in_one_folder"
setwd(wd)
Then you get the list of files you want from your folder
files <- list.files()
Put into a list with a loop
dflist <- list()
for(i in 1:length(files)){
 dflist[[i]] <- read.csv(files[i])
}

Then use tidyverse solution
library(tidyverse)
map(dflist, ~.$GeneID) %>% #creates a list just of the column of interest
    reduce(intersect)  #applies the intersect function cumulatively to the list

If you need a package it would be better if you state which ones. Also with readr package you can read already a bunch of csv files in one single call

I mentioned tidyverse, but I just added library(tidyverse) in the example to be extra clear. Also, I felt as though a loop example would be better for explaining how to construct them..

" a loop example would be better for explaining how to construct them" this is what im looking for even though i dont use much loop etc because Im sure or my concepts are pretty unclear about it.

@h3ab74  i ran your code    map(dflist, ~.$GeneID) %>% reduce(intersect)
I get this NULL .Am i doing something wrong

the part of the code with ~.$GeneID is supposed to be the column name of all your dataframes in the list. So, if the column with the gene names in your dataframes does not have the name 'GeneID', change the GeneID part in the code to one that works with your dataframes

In my case the column name of all my data frame header is Gene, so i ran this
map(dflist, ~.$Gene) %>%
reduce(intersect)
character(0)  I get character(0)

Can you edit your question to include an example of how the dataframes are structured so I can see why you're getting nothing? Also, please note that this solution will get you the genes that are similar in all of the dataframes

@h3ab74 these are my files https://drive.google.com/open?id=1laieKec1sr4rWa4SIhUvHcLSerrWf4OC my data frame which i import is single column with a header name Gene

it worked ,it had name issues.


text:
K mean clustering issue
So I have expression data set with 4 healthy sample and 4 disease, to determine the number of K I used mclust which i ran on the list of differentially expressed genes so I got around 7 cluster from Mclust which im using as an input for my K-mean .
Now my data set is like all the differentially expressed genes it contains both UP and down regulated genes so basically it may be UP in my healthy and Down in patient and vice-versa ,4 likely possibilities.
So now I have this doubt that from that I'm getting 7 cluster.Is that result is correct ?

What do you mean that you have 4 possibilites? Which constrasts did you do to calculate the differentially expressed genes? Did you compare healthy to disease? Then you only have 2 options.

There is generally no such thing as correct or incorrect when it comes to clustering. When you specify $K$, the program did it's statistical algorithms. Unless the the program had a bug, your results were statistically valid.
You will need to check the clustered labels visually. Only your biological knowledge can tell you the optimal $K$.

7 clusters may be correct in that it's the optimal number according to the Bayesian information criterion (or whatever value you had mclust optimize), but that doesn't make it necessarily the most sensible value. Have a look at your data and visually look at what might make sense. It's likely that the clusters are subdividing up and down regulated genes into various subgroups, which may or may not be useful, depending on your goals.

yes I used BIC ,and this is what my PI was asking if those cluster are even correct or make sense ...may be I would try some other clustering and see if get similar number I would let know...

As you're learning, clustering tends to be an iterative process.

very naive question when i define the number of k , the k mean clustering happens on my rows or columns or both ?

Rows are clustered using the values in the columns.


text:
Compare nucleotide mutation with amino acid
To begin with I have two nucleotide sequence to compare I want to see how they might have evolved in terms of sequence similarity and differences. One of the way i can think of is align the sequences. But i m not getting how do i associate it with mutation.
Now question is how do i know about mutation. I have the respective amino acid data from which i know there been mutation which is making changes how the binding happens between two sequences.
So now i have to associate nucleotide mutation with amino acid sequence data.
Any suggestion or help would be really appreciated. As always if there is a R way to performing the analysis then it would be really good.

Are you familiar with non-synonymous vs synonymous analysis and is this part of the phangorn package?

i did come to read that let me explore it..will get back

Thats probably what you want ... depends on your investigation

so my simple idea is compare the respective nucleotide sequence and see if there exist a mutation if yes then do the same in the protein sequence. So to see mutation I'm not sure how to start ..alignment perhaps not sure though..how to approach it..

You must have done an alignment to generate the phangorn tree and if you have branches and stuff you must have mutations. An alignment viewer that converts in and out of amino acid translation is a place to start

yes i have done mutation and i have the tree as well ."An alignment viewer that converts in and out of amino acid translation" jalview would be a good start i guess..

Cool, I will sketch a few ideas below when I have a minute

that would be really kind of you

I  would output the aligned nucleotide sequences and input them into an alignment editor such as Sequtron (Mac), MEGAX(particularly PC) or Jalview (Linux).
The alignment editors will allow you to flick between the triplet codon and the amino acid sequence. Therefore spotting which triplet codon mutations associate with each given amino acid mutation is easy.There is no R equivalent that I am aware of.
More generally I would do a non-synonymous / synonymous analysis which fangorn will perform. This assesses the number of mutations in which cause an amino acid change against the number of mutations that cause a 'silent' mutation (no amino acid change). A synonymous phylogeny will represent all mutations which don't change the tree. A non-synonymous phylogeny is the same as an amino tree.
There are character mapping approaches which will map an amino acid mutation onto a phylogeny as an alternative.

I did run the kaks but ran into this https://bioinformatics.stackexchange.com/questions/13499/calculating-synonymous-and-non-synonymous

Hi @krushnachChandra I've answered your question. Why it fails is easy, but correcting it is more complication, however this is needed before any form of amino acid to nucleotide mutation analysis can take place.


text:
Finding out-group in unrooted tree
Im trying to find the out-group in the alignment my data files file a and  file b both the case i'm getting NA when i run this     pars.a$tip.label[pars.a$edge[which.max(pars.a$edge.length),2]]
I tried NJ instead of upgma i do still get NA not sure what is the issue ,but at the same time when i used other alignment file with 50 sequence  works fine.
a<-read.alignment("abhi_seq/final_abhi_tempa1.fas", format="fasta")
a.phydat<-as.phyDat(a)
dist.a.phydat<-dist.dna(as.DNAbin(a.phydat))
upgma.a<-upgma(dist.a.phydat)
parsimony(upgma.a,a.phydat)
pars.a <- optim.parsimony(upgma.a, a.phydat)
pars.a<-acctran(pars.a, a.phydat)
pars.a$tip.label[pars.a$edge[which.max(pars.a$edge.length),2]]

the fact is that the present upgma tree as a midpoint rooting here, corresponding to node 46.
getRoot(upgma.a)

so for this kind of tree, I recommend you to use the midpoint() rooting function.
pars.a <- optim.parsimony(upgma.a, a.phydat)
pars.a<-acctran(pars.a, a.phydat)
plot(midpoint(pars.a))

let me try it out and how do we decide which rooting to use? based on observation after plotting

yes visual inspection seems relevant. It totally depends of your input data, whether your sequences are strongly conserved (without evident outgroup) or more divergent. I think you may write a control loop after constructing the upgma tree, by checking if the ouptput of getRoot(upgma.a) returns a node or a tip. If it is a node, then it might proceed with the midpoint() version.

I haven't tried with ML but this mid point should work with ML as well?

yes, probably! it is probably a more universal solution, but not the necessarily the best one!

last time i was looking to compare tress which you had suggested turn them into dedrogram and then for comparison you suggested something this '
pars.a.rooted<-root(pars.a, outgroup="ABN10911", resolve.root=T)
pars.a.rooted.dd<-as.dendrogram(force.ultrametric(pars.a.rooted))' how to do the same  here shall i put node 46 as the outgroup?

It won't work because you have only 45 tip labels, if I remember well. I think you may run: pars.a.rooted.dd<-as.dendrogram(force.ultrametric(midpoint(pars.a)))

works perfectly fine for parsimony and...for pars for ML i m doing this rectify me if im wrong "ml.a.rooted.dd<-as.dendrogram(force.ultrametric(midpoint(fit$tree)))
"

seems ok. It is true that in the previous code, we used  root(midpoint(fit$tree), outgroup="AAA64460", resolve.root=T) to force using a given tip as root. Using force.ultrametric(midpoint(pars.a)) should work too indeed.

Let us continue this discussion in chat.


text:
Chromatin accessibility level on promoters with different CpG ratio calculation
I want to categorize promoters based upon high, intermediate and low-CpG content (high-CpG-density promoters (HCPs), intermediate-CpG-density promoters(ICPs), and low-CpG-density promoters (LCPs)).
So the data I have is for promoter which i have annotated 1000+/- around TSS and taken them as something less than 1kb is promoter region and beyond distal.
In term of tool i have used chipseeker to annotate and do the above step.
So now if i have to find the CpG density as i have mentioned above how do i get that?

you can either calculate the GC content? And you can also check overlap with annotated cpg island from UCSC... Not very sure what is the question

can i get the GC content from  atac seq data ? if yes..how that is my question I mean what input do i require

did i get you wrong? It seems like you have annotated peaks. It's a matter of getting the sequences in the annotated peaks and calculating the number of CGs (sorry I typed too quickly previously, it should be CGs instead of GC content)

"It seems like you have annotated peaks." yes correct. "It's a matter of getting the sequences in the annotated peaks and calculating the number of CGs" okay that looks straight forward though can it be done in R ? can you suggest me some library

Below I import the peaks call in narrowbed format from macs2 into a GRanges object. If you have a bed file, you can just use rtracklayer for it:
library(rtracklayer)
library(BSgenome.Hsapiens.UCSC.hg38)

peaks_narrowbed = "https://www.encodeproject.org/files/ENCFF846JAO/@@download/ENCFF846JAO.bed.gz"

extraCols <- c(signalValue = "numeric", pValue = "numeric",
                          qValue = "numeric", peak = "integer")
gr <- import(peaks_narrowbed, format = "BED",extraCols = extraCols)

head(gr)

GRanges object with 6 ranges and 6 metadata columns:
      seqnames        ranges strand |
         <Rle>     <IRanges>  <Rle> |
  [1]     chr1   10049-10326      * |
  [2]     chr1 180686-181040      * |
  [3]     chr1 186648-186903      * |
  [4]     chr1 191293-191604      * |
  [5]     chr1 267872-268146      * |
  [6]     chr1 586110-586324      * |
                                                   name     score signalValue
                                            <character> <numeric>   <numeric>
  [1] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_1       135     7.10907
  [2] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_2       309    11.65888
  [3] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_3        81     5.40289
  [4] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_4        56     4.54981
  [5] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_5       145     7.39343
  [6] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_6       116     6.54035

Then we extract the sequences that are within these regions, you can of course redefine the region by extending it from the peak:
library(Biostrings)
SEQ = getSeq(Hsapiens,gr)
head(SEQ)

A DNAStringSet instance of length 6
    width seq
[1]   278 TAACCCTAACCCTAACCCTAACCCTAACCCTAAC...CCCCAACCCCAACCCCAACCCCAACCCCAACCC
[2]   355 CTAACCCCTAACCCCTAACCCTAACCCTACCCTA...CTCAGCCGGCCCGCCCGCCCGGGTCTGACCTGA
[3]   256 GCAACCACCTGAGCGCGGGCATCCTGTGTGCAGA...CCCCTTCTTTCCATTGGTTTAATTAGGAACGGG
[4]   312 CAATCAGCAGGGACCGTGCACTCTCTTGGAGCCA...CCCTGCCCCTGTCTCTTCCGTGCAGGAGGAGCA
[5]   275 TGGATGTGTGCATTTTCCTGAGAGGAAAGCTTTC...TTGTTTAGTTTGCGTTGTGTTTCTCCAACTTTG
[6]   215 GAGGCTGAAGGAGACTGATGTGGTTTCTCCTCAG...CCAACTTTGTGCCTCATCAGGAAAAGCTTTGGA

This is the CG, and you can decide how you want to bin your regions:
freq = dinucleotideFrequency(SEQ,as.prob=TRUE)
hist(freq[,"CG"],br=100)

thank you veryyyyy much..i did look into couple of papers but it didn't occur to me how they did it..i will try this and update you..

"This is the CG, and you can decide how you want to bin your regions:" very naive  question the sequence also contain the "AT " how do you define CG?

sorry I typed too quickly, so in the code above, the result from dinucleotideFrequency() contains the proportion of dinucleotides that are CG, which is what you need, if you are interested in CpG?

yes i need CG for example this "TGGATGTGTGCATTTTCCTGAGAGGAAAGCTTTC"  "the result from dinucleotideFrequency() contains the proportion of dinucleotides that are CG" this sequence which i have put also contains AT my question is is the function "dinucleotideFrequency" only look for CG only??

of course not.. see when i wrote this freq = dinucleotideFrequency(SEQ,as.prob=TRUE), freq is a matrix. the next line i called out freq[,"CG"], did you look at the matrix?

" freq[,"CG"]" im extremely sorry not sure how did i miss that!!!!  now i got it and big thank to make it clear to me...


text:
Comparing aligned amino acids to codon
I have a set of 4 amino acids which i have aligned,but i want to compare then with respective  nucleotide sequence in terms of change codon level.
I have the alignment file where i want to take only for the C-terminal domain region region only the CTD starts with amino acids “NITNLC” till it ends with “HAPATV” and rest of the sequence i dont want to take.
In short i want to take the part sequence which starts from "NITNLC" to "HAPATV" then i would like to compare what are the changes in the nucleotide sequence at codon level.

To start with I have used the seqinr library to generate nucleotide alignment from aligned protein sequence.
reverse.align {seqinr}

Then the next step that i want to compare each of my amino acid sequence to my reference sequence which in my case is P0DTC2to rest of the sequence for example I want to compare P0DTC2 to K9N5Q8to the above mentioned CTD region which starts from NITNLC and ends with HAPATVand find what are the codon level changes if there is any I would like to report both in amino acid level and codon level.

First part I guess my approach is correct. Second part I'm not sure how to proceed i guess it's much more than simple parsing I guess!!.
Any help or suggestion would be really appreciated ,and if there exist a R based solution that would be really welcome.
Files
Amino acid file nucleotide file Aligned amino acid file Reverse alignment nucleotide file using seqinr

To set the ball rolling, I tried to implement it manually with your files. Overall, there are three steps:

Locate NITNLC (or HAPATV) in protein P0DTC2.
Locate mismatched amino acids between proteins P0DTC2 and K9N5Q8 within the range from step 1.
Print the amino acid and the DNA codon of both proteins from step 2.

It works, but only for the first 60 amino acids. I wonder if it's due to the amino sequences in AMINOOO_seq_removed.fasta repeat itself every 60 acids. But why?
#For example, the first three lines of protein P0DTC2
>P0DTC2
MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFS MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFS
NVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIV NVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIV
NNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLE NNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLE
...

Step 0: Read the files.
library('seqinr')

#align file containing protein sequences
count_added <- read.alignment('count_added_.clustal_num', format='clustal')
names(count_added$seq) <- count_added$nam

#DNA sequences
rev3.aln <- read.alignment('rev3.aln', format='fasta')
names(rev3.aln$seq) <- rev3.aln$nam

Step 1: Locate NITNLC (or HAPATV) in protein P0DTC2. Again, there are two repeats of NITNLC 60 acids apart (at 807 and 870).
locate <- function(seq, find)
  {address <- gregexpr(paste(strsplit(find, '')[[1]], collapse='[^a-z]*'), seq)
  #substr(seq, address[[1]][1], address[[1]][1]+attr(address[[1]], 'match.length')[1]-1)

  return(list(start=as.numeric(address[[1]]),
             end=as.numeric(address[[1]] + attr(address[[1]], 'match.length') - 1)))
  }

locate(seq=count_added$seq[['P0DTC2']], find='nitnlc')
#start
#807 870
#end
#812 875
locate(seq=count_added$seq[['P0DTC2']], find='hapatv')
#start
#1212 1279
#end
#1217 1285

Step 2: Locate mismatched amino acids between proteins P0DTC2 and K9N5Q8 within the range from step 1. Instead of using the range 812 - 1279, I chose range 1 - 20 for demonstration purposes.
compare <- function(seq1, seq2, after=0, before=100000)
  {seq1_ = strsplit(seq1, '')[[1]]
  seq2_ = strsplit(seq2, '')[[1]]
  ind = which(seq1_ != seq2_ & grepl('[a-z]',seq1_) & grepl('[a-z]',seq2_))
  ind = ind[ind>after & ind<before]
  #seq1_[ind[1]]
  #seq2_[ind[1]]

  return(ind)
  }
compare(seq1=count_added$seq[['P0DTC2']], seq2=count_added$seq[['K9N5Q8']], after=1, before=20)
# [1] 5  7  8  9 10 13 14

#protein comparison
#K9N5Q8      MIHSVFLLMFLLTPTESYVD
#P0DTC2      ----MFVFLVLLPL------
#<mismatch>      5 7890  34

Step 3: Print the amino acid and the DNA codon of both proteins from step 2. Note that ind is based on the align file.
print_amino_codon <- function(ind, seq, seq_gene)
  {locate_amino <- gregexpr('[a-z]', seq)[[1]]
  if (!ind %in% locate_amino) return(NA)
  ind2 = match(ind, locate_amino)

  return(c(amino=substr(seq, ind, ind), codon=substr(seq_gene, ind2*3-2, ind2*3)))
  }

codon(ind=5, seq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
#amino codon
#  "v" "atg"
codon(ind=6, vseq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
#amino codon
#  "f" "ttc"
codon(ind=7, seq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
#amino codon
#  "l" "ttg"

codon(ind=5, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
#amino codon
#  "m" "atg"
codon(ind=6, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
#amino codon
#  "f" "ttt"
codon(ind=7, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
#amino codon
#  "v" "gtt"

#protein comparison
#K9N5Q8      MIHSVFLLMFLLTPTESYVD
#P0DTC2      ----MFVFLVLLPL------
#<mismatch>      5 7890  34
#<print>         ^^^

#K9N5Q8 gene codon
#gtg ttt cta ctg atg ttc ttg tta aca
#                ^^^ ^^^ ^^^
#P0DTC2 gene codon
#atg ttt gtt ttt ctt
#^^^ ^^^ ^^^

Perfect @RyanSYKwan its alot of code

thanks a lot this will be really really helpful for my future work

Glad it helps! I learn a lot from you too.


text:
Chip-seq down-sampling problem
Control sample chip read stat
56234445 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
55173337 + 0 mapped (98.11% : N/A)
56234445 + 0 paired in sequencing
28078024 + 0 read1
28156421 + 0 read2
53808040 + 0 properly paired (95.69% : N/A)
54626508 + 0 with itself and mate mapped
546829 + 0 singletons (0.97% : N/A)
162718 + 0 with mate mapped to a different chr
112602 + 0 with mate mapped to a different chr (mapQ>=5)

IgG sample chip seq read stat
32745255 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
31546173 + 0 mapped (96.34% : N/A)
32745255 + 0 paired in sequencing
16362391 + 0 read1
16382864 + 0 read2
30934194 + 0 properly paired (94.47% : N/A)
31269484 + 0 with itself and mate mapped
276689 + 0 singletons (0.84% : N/A)
79786 + 0 with mate mapped to a different chr
49658 + 0 with mate mapped to a different chr (mapQ>=5)

knockout sample Rep 1 stat
28044731 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
27459222 + 0 mapped (97.91% : N/A)
28044731 + 0 paired in sequencing
14014483 + 0 read1
14030248 + 0 read2
26685254 + 0 properly paired (95.15% : N/A)
27162846 + 0 with itself and mate mapped
296376 + 0 singletons (1.06% : N/A)
116684 + 0 with mate mapped to a different chr
88965 + 0 with mate mapped to a different chr (mapQ>=5)

Knockout sample rep 2 stat
32410094 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
31774629 + 0 mapped (98.04% : N/A)
32410094 + 0 paired in sequencing
16196814 + 0 read1
16213280 + 0 read2
30802808 + 0 properly paired (95.04% : N/A)
31561218 + 0 with itself and mate mapped
213411 + 0 singletons (0.66% : N/A)
89688 + 0 with mate mapped to a different chr
63838 + 0 with mate mapped to a different chr (mapQ>=5)

Do i need to down-sample with respect to the IgG for control sample which has lot more reads at the same time Knockout sample 1 has less read. So is it good to go for further analysis or i need to down sample them ?
Any suggestion is really appreciated

Usually downstream programs will take sequencing depth differences into account for you, so you only need to downsample if you get really funky results.

"if you get really funky results." i did plot using your deeptools i do see a small peak in IgG  so i feel there might be some issue.. which one would be better IgG or input im not clear .for the comparison ..any insight if IgG would be better or input ?

control is almost always better than IgG. For this reason don’t bother with IgG.


text:
Pausing index calculation
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))

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

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

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

it's better to spend time on something more interesting like what's in the promoter or SUPER enhancer, than thinking about pausing index

"Warning message:
Computation failed in stat_signif():
not enough 'y' observations " i did get this as im using full data frame

damn i just copy pasted your code and fail to see that...............

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?

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


text:
Module preservation and hub genes finding
I ran module preservation analysis between two condition and on the basis of preservation z summary i would to go for downstream analysis.So I m interested in looking for hub-genes of those modules which have lowest z summary. But so far in the tutorial i didn't find a section in the WGCNA tutorial where how to find hubgenes after module preservation analysis was mentioned.
Any suggestion would be really appreciated

Simply select the module(s) you are interested in and look for hub genes using the standard calculations. Having done a module preservation analysis does not change the procedure in the slightest.

but so far what i used to do i run network analysis on individual condition and find hub genes but in module preservation i haven;t found the step or the output from which I can look for hub-genes it would be really helpful ,if you can cite me "using the standard calculations. " in your site

Well, if you ran module preservation, you must have already previoulsy done a network analysis of at least the reference data. Go back to that network analysis and focus on the modules flagged by the preservation analysis. Use kME in those modules to select the top hubs. Section 3d of Tutorial I shows you how to create a table of kMEs.

thank you for the reference


text:
Peak-calling using homer
I have a total of 78 ATAC seq samples from which Im trying to do a peak call .I tried this batchParallel.pl findPeaks peaks function but i couldn't find the output what it makes.
So I tried using loop but again the peak file is not written due to incorrect argument usage I guess.
These are my sample directories where i have all the tags
SRR2920564  SRR2920565

I used homer to make tag directories and now i using it for peakcalling.
here is my command.
for d in */ ; do findPeaks $d -style factor -o "${d}.txt" ; done

What i need is all the peakfile should be written as the directory name.
Such as if directory is "SRR2920564" then my peak file should be written in the same directory as SRR2920564.txt
I do run the loop it runs successfully but no output written. May be a mundane mistake am not able to figure it out.
Any suggestion or help would be highly appreciated .

Your output file is a bit weird.
findPeaks directoryname/ -style factor -o directoryname/.txt
You would be attempting to create a file called ".txt" inside each directory? (Also do you need those quotes).
Did you try running this with "-o auto" and checking in each directory if the corresponding files (i.e peaks.txt, regions.txt, transcripts.txt, etc..) were generated?

I tried running auto it gives an output peak.txt but when it would be easy if i can name them as the directory name for as unique name for each file

You just need to strip the / from $d, for which there are a number of options:
for d in */ ; do
    dname=`basename $d`
    findPeaks $d -style factor -o ${dname}.txt
done

or
for d in */ ; do
    dname=${d%%/}
    findPeaks $d -style factor -o ${dname}.txt
done

There are likely other ways one could go about this. Having said that, since you're not running these in parallel you could simply move the default output file (peak.txt) to a file name of your choice (using one of the strategies above).


text:
Tree cut issue in WGCNA
library(WGCNA);
# The following setting is important, do not omit.
library(flashClust)
enableWGCNAThreads()
options(stringsAsFactors = FALSE);
#Read in the female liver data set
femData = read.csv("CTRL_WGCNA.txt",sep = "\t");

head(femData)
# Read in the male liver data set
maleData = read.csv("BPD_WGCNA.txt",sep = "\t");
# Take a quick look at what is in the data sets (caution, longish output):
dim(femData)
names(femData)
dim(maleData)
names(maleData)

#=====================================================================================
#
#  Code chunk 2
#
#=====================================================================================

# We work with two sets:
nSets = 2;
# For easier labeling of plots, create a vector holding descriptive names of the two sets.
setLabels = c("Female liver", "Male liver")
shortLabels = c("Female", "Male")
# Form multi-set expression data: columns starting from 9 contain actual expression data.
multiExpr = vector(mode = "list", length = nSets)

multiExpr[[1]] = list(data = as.data.frame(t(femData[-c(1)])));
names(multiExpr[[1]]$data) = femData$Gene;
rownames(multiExpr[[1]]$data) = names(femData)[-c(1)];
multiExpr[[2]] = list(data = as.data.frame(t(maleData[-c(1)])));
names(multiExpr[[2]]$data) = maleData$Gene;
rownames(multiExpr[[2]]$data) = names(maleData)[-c(1)];
# Check that the data has the correct format for many functions operating on multiple sets:
exprSize = checkSets(multiExpr)

#=====================================================================================
#
#  Code chunk 3
#
#=====================================================================================

# Check that all genes and samples have sufficiently low numbers of missing values.
gsg = goodSamplesGenesMS(multiExpr, verbose = 3);
gsg$allOK

#=====================================================================================
#
#  Code chunk 4
#
#=====================================================================================

if (!gsg$allOK)
{
  # Print information about the removed genes:
  if (sum(!gsg$goodGenes) > 0)
    printFlush(paste("Removing genes:", paste(names(multiExpr[[1]]$data)[!gsg$goodGenes],
                                              collapse = ", ")))
  for (set in 1:exprSize$nSets)
  {
    if (sum(!gsg$goodSamples[[set]]))
      printFlush(paste("In set", setLabels[set], "removing samples",
                       paste(rownames(multiExpr[[set]]$data)[!gsg$goodSamples[[set]]], collapse = ", ")))
    # Remove the offending genes and samples
    multiExpr[[set]]$data = multiExpr[[set]]$data[gsg$goodSamples[[set]], gsg$goodGenes];
  }
  # Update exprSize
  exprSize = checkSets(multiExpr)
}

#=====================================================================================
#
#  Code chunk 5
#
#=====================================================================================

sampleTrees = list()
for (set in 1:nSets)
{
  sampleTrees[[set]] = hclust(dist(multiExpr[[set]]$data), method = "average")
}

#=====================================================================================
#
#  Code chunk 6
#
#=====================================================================================

pdf(file = "SampleClustering.pdf", width = 12, height = 12);
par(mfrow=c(2,1))
par(mar = c(0, 4, 2, 0))
for (set in 1:nSets)
  plot(sampleTrees[[set]], main = paste("Sample clustering on all genes in", setLabels[set]),
       xlab="", sub="", cex = 0.7);
dev.off();

  #=====================================================================================
#
#  Code chunk 7
#
#=====================================================================================

# Choose the "base" cut height for the female data set
baseHeight = 80
# Adjust the cut height for the male data set for the number of samples
cutHeights = c(80, 100*exprSize$nSamples[2]/exprSize$nSamples[1]);
# Re-plot the dendrograms including the cut lines
pdf(file = "SampleClustering.pdf", width = 12, height = 12);
par(mfrow=c(2,1))
par(mar = c(0, 4, 2, 0))
for (set in 1:nSets)
{
  plot(sampleTrees[[set]], main = paste("Sample clustering on all genes in", setLabels[set]),
       xlab="", sub="", cex = 0.7);
  abline(h=cutHeights[set], col = "red");
}
dev.off();

  #=====================================================================================
#
#  Code chunk 8
#
#=====================================================================================

for (set in 1:nSets)
{
  # Find clusters cut by the line
  labels = cutreeStatic(sampleTrees[[set]], cutHeight = cutHeights[set])
  # Keep the largest one (labeled by the number 1)
  keep = (labels==1)
  multiExpr[[set]]$data = multiExpr[[set]]$data[keep, ]
}
collectGarbage();
 # Check the size of the leftover data
exprSize = checkSets(multiExpr)
exprSize

$nSets
[1] 2

$nGenes
[1] 8213

$nSamples
[1] 0 0

$structureOK
[1] TRUE

I think the issue is with the tree cut height because when i check the label all of them are zeroes i guess that is issue so far which is why im getting exprSize nsamples as 0 0instead of my actual sample size.
I have tried it with different cut height but so far no success may its a pretty simple issue but im not able to figure it out.
Data which im using here is the link
Any suggestion or help would be highly appreciated

It is likely that the cut heights set in the code (80 and 100*ratio of sample numbers) are too low. Look at the sample trees plotted into a pdf file (SampleClustering.pdf) and note the merging heights. You need to set the cut heights appropriately so that most samples go into a single cluster.
Peter

"(SampleClustering.pdf) and note the merging heights. You need to set the cut heights appropriately so that most samples go into a single cluster." so my y axix upper limit is 120 .Should i set it to 120 as its the highest height .I will attach the fig which was generate for the clustering

The issue is resolved now


text:
WGCNA module preservation analysis
I have case control studies,4 different condition where I have 1 control and against which I'm have 3 comparisons to make.
I did WGCNA for individuals samples .Here I have the kept the number of genes equal in all the cases .
Now my question is whether the number of modules which are obtained from WGCNA would be similar in number if I run module preservation between two condition .
Lets say i get 15 modules in Control condition and 16 treated condition if I run module preservation the number of modules which would come would be 15+16= 31 or its independent .
The power I have kept same for both WGCNA and module preservation analysis. But Still I get different module numbers.
Any suggestion would be highly appreciated

Your question is not clear, you seem to mix the relevant terminology...
First, you cannot run WGCNA for individual samples, only for individual conditions (with or without WT, depending on what question you're trying to answer and assuming you have enough samples in each condition). If that's what you have done, fine. The number of modules in each independent analysis will generally be different, that's because each condition is different. The module identification built into WGCNA as the default approach does not allow the user to specify the number of modules beforehand, but if you insist on getting a pre-specified number of modules, you could run WGCNA using the "manual" way (for example, Tutorial I, section 2b) and substitute another clustering method (e.g., PAM) that allows you to specify the number of clusters. You may find the results inferior though.
Second, module preservation does not result in new modules. Module preservation evaluates whether modules found in data set A are preserved in a different data set B. To get new, common modules, you could run a consensus network analysis on pairs of conditions or all 3 conditions together. The number of modules you would find in a consensus analysis cannot be predicted from the number of modules you find in individual data sets.

"assuming you have enough samples in each condition)." sorry yes i mixed up and yes i have control have like 55 samples and treatment like 50 samples .Thank you for the clarification now it seems bit more clear I will read the tutorial a bit more

Conceptual clarification I have done as I said Control sample WGCNA as well as treatment ones with numerous biological replicate .I get various modules and yes "umber of modules in each independent analysis will generally be different, that's because each condition is different." yes this is what i get .So if do module preservation analysis where I want to see if my "data set A are preserved in a different data set B" then the number of resulted module would be same as what I did WGCNA for Control sample where I would see w.r.t to my reference. Am I explaining it correctly ?

Sorry, I still don't understand your question. Module preservation outputs, for each module in a reference network (A), preservation statistics of the module in a test network (B). Module preservation does not change the number of modules. Some modules from A may be preserved in B and some not; sometimes all modules are preserved, sometimes none of the modules are preserved. The number of modules from A that are preserved in B may be smaller or larger than the number of modules found in a separate network analysis of B.

Thank you for the clarification


text:
bedtool intersect multiple bed files
I have 7 biological replicate of a normal sample after peak-calling I converted it into bed files and now I want to find what are the common region between those .
So i tried with bedtool intersect
my command is this
bedtools intersect -a SRR2920506.filtered.bed -b SRR2920531.filtered.bed SRR2920466.filtered.bed SRR2920478.filtered.bed SRR2920507.filtered.bed SRR2920532.filtered.bed  > answer.bed

it says " the intersect tool can detect overlaps between a single -a file and multiple -b files (instead of just one previously). One simply provides multiple -b files on the command line."
My doubt is the file that goes to "a" argument is not query which i have to detect overlaps against "b" files as these are all biological replicate. Am i doing the right way?
Secondly I tried with bedops but the number of overlaps which is reported is really low when i compare the result which i get from bedtools intersect function.
Any suggestion or help regarding how to intersect multiple would be really appreciated

If you want to find the intersection between ALL your bedfiles, you can try multiIntersectBed (available since bedtools 2.14.3).
It should work like this :
multiIntersectBed -i SRR2920506.filtered.bed SRR2920531.filtered.bed SRR2920466.filtered.bed SRR2920478.filtered.bed SRR2920507.filtered.bed SRR2920532.filtered.bed  > answer.bed

To sort your bedfiles by chromosome then by start :
sort -k 1,1 -k2,2n myfile.bed

ERROR: input file: (SRR2920466.filtered.bed) is not sorted by chrom then start. this came as i tried with multiIntersectBed

Yep this is described as the second line of the usage. I just updated the answer


text:
multiBamSummary output chromosome region annotate issue
I was looking how to take readcount from ATAC seq data , this post which suggested to use bedtools and featurecounts to take out readcounts .But multiBamSummary seems to be pretty straight forward. Now I have my count file but only issue is there are no gene names. So now is it possible to annotate these chromosome regions? if so how ?
Any suggestion or help would be highly appreciated

Assuming you used a BED file as input then the output is in the same order, so just merge the two files together.

i used bam files in case of bed files but yes the orders are same .But would that help me in annotating the region ?

I assume you used multiBamSummary BED-file, did you use multiBamSummary bins instead?

"multiBamSummary bins" yes that is what i used what i saw in the manual

They're bins of the size you specified that tile the genome then.

okay.Now any way to annotate those regions ? the intervals which i see in the output ,as i would like to take that for downstream analysis where i would like to see the gene name

Overlap/intersect them with a GTF file.

Conceptual doubt ,would it be similar to annotating the peaks like that of homer annnotate peak function ? or I can simply map the the chromosome start end to gtf file.

"Conceptual question" (a doubt is generally not something one "has" in English). I'd just simply map things to the GTF file.

thank you @Devon clarifying my question and the English


text:
DNA bioinformatic search and methylation analysis by qMSP
Im quoting the methods which was used

We performed a bioinformatics analysis to identify potential CpG
  islands in the genomic sequences of RANKL and OPG genes. The USCS
  Human Genome Browser public database (genome.ucsc.edu), the Primer
  Express software (available from Applied Biosystems), and CpG Island
  Explorer software (cpgie.sourceforge.net) were used to explore the
  nucleo- tide sequences surrounding de transcription start site (TSS)
  of RANKL isoforms I and II (NG_008990.1 Ref Seq Gene), and OPG
  (NG_012202.1 Ref Seq Gene). A (C + G)/total bases ratio

0.5 and CpG observed/CpG expected ratio >0.6 were used as criteria.

Im looking for the region which they have used
This is the paper
"Role of DNA methylation in the regulation of the
RANKL-OPG system in human bone"

I tried using ucsc to check the region but didn't get much. The idea is to find the region to design the primer for the same
How can i get the region from ucsc ,the CpG island around the TSS of the above mentioned genes which are RANKL and OPG
Any suggestion or help would be highly appreciated

It's likely that they accomplished this by using the UCSC Table Browser to:

Find the loci of the two isoforms you mention
Finding CpG islands (Group = Regulation, Track = CpG Islands) near those transcription start sites
Design primers using Applied Biosystem's software

Alternatively, they list the exact primers used in Table 1.
They don't detail exactly how they did it, so if you have specific questions it might be worth it to email the corresponding author

thank you for the suggestion i would try it out and yes it would be better if i drop them a mail


text:
samtools fixmate threads for uncompressed bam data in pipeline
Samtools' documentation regarding duplicate marking suggests the following pipeline:
samtools collate -@ 4 -O -u example.bam | samtools fixmate -@ 4 -m -u - - | samtools sort -@ 4 -u - | samtools markdup -@ 4 - markdup.bam

With the following comment:

The -u option produces uncompressed output. Since the output is not being written to disk then size is unimportant and avoiding compress/decompress reduces the amount of work that needs doing. The -@ option turns on threading which can speed things up.

However, the documentation of samtools fixmate (*) says:

-@, --threads INT
Number of input/output compression threads to use in addition to main thread [0].

Since option -u is used upstream of fixmate and since fixmate also uses option -u, I expect fixmate to have no "input/output compression" to perform at all. Is providing threads to fixmate going to be of any help under these circumstances?
(*) http://www.htslib.org/doc/samtools-fixmate.html
(samtools collate will read compressed input, samtools sort also uses threads for sorting, and samtools markdup will compress its output, therefore I expect some effect for option -@, --threads there.)

I think allowing samtools fixmate to use additional threads is unlikely to provide any performance benefit given that -u (for uncompressed output) eliminates the need for compression/decompression. I suspect that the main computational effort for samtools fixmate (i.e. fixing mate-pair information) will remain single-threaded. That said, the best way to verify any and all performance assumptions is to test and profile the pipeline. In my experience though, specifying -@ (for additional threads) is most beneficial when used with samtools sort where multi-threading can significantly improve performance. You may not even notice the impact of additional decompression or compression threads at the beginning and end of the pipeline.


text:
Adding entries to bigwig file
I generate bigwig files using a shell script based on bedtools genomecov (to generate a bedgraph from a bam file) bedmap (to compute means across 10 bp bins) and bedGraphToBigWig to convert the binned bedgraph into bigwig.
Sometimes, the bam file has no data on some chromosomes, this results in a bedgraph file with no entries for the empty chromosomes.
I would like to add the missing entries. How could I do this?
I tried to see what I could do with pyBigWig, opening the bedgraph file in "a" mode, but don't really know how to proceed:
# I have the chromosome sizes for my genome of interest already in a list:
chrom_sizes
#[('I', 15072434),
# ('II', 15279421),
# ('III', 13783801),
# ('IV', 17493829),
# ('MtDNA', 13794),
# ('V', 20924180),
# ('X', 17718942)]
bw_test = pyBigWig.open("path/to/bigwig", "a")
bw_test.chroms()
# {'IV': 17493829, 'V': 20924180, 'X': 17718942}
bw_test.addHeader(chrom_sizes)
# RuntimeError: Received an error in bwCreateHdr

Presuming you still have your BAM file sitting around somewhere:
#!/usr/bin/env python
import pysam
import pyBigWig

bam = pysam.AlignmentFile("alignments.bam")
ifile = pyBigWig.open("your current file.bw")
ofile = pyBigWig.open("your new file.bw", "w")

# Create/Add the header
hdr = []
for chrom in bam.header['SQ']:
    hdr.append((chrom['SN'], chrom['LN']))
bam.close()
ofile.addHeader(hdr)

# Write the bigWig
for _ in hdr:
    try:
        ints = ifile.intervals(_[0])
        starts = [x[0] for x in ints]
        ends = [x[1] for x in ints]
        vals = [x[2] for x in ints]
        ofile.addEntries([_[0]] * len(starts), starts, ends=ends, values=vals)
    except:
        # intervals() will throw an error if the chromosome isn't present
        pass
ofile.close()
ifile.close()

The reason you can't append to a bigWig file is that it has a various indices and zoom levels typically built in. You'd have to add to them without overwriting neighboring blocks of data. In theory this might sometimes be theoretically possible, but not typically and I'm not aware of tools that support this.

If I understand correctly, the bam file is not necessary in my case, because I already have the chrom_sizes list, which should be the same as the hdr you obtain by parsing the bam header.

Correct, if you already have such a list then you don’t need a BAM file.

Just a detail: It seems strange to me to use a _ placeholder as looping variable if you then use it to access its first element. I suppose for (chrom, _) in hdr: would work as well while avoiding you the _[0].

Yup that’d work too.

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
How can I easily get the read size distribution of reads mapping on a certain set of regions?
Suppose I have a BAM file indicating where reads in a library have mapped, and a bed file describing a set of genomic regions.
Is there a way to easily get the size distribution of the reads mapping on this set of regions?

This can be done rather simply by combining this samtools based answer and bioawk:
samtools view -L regions.bed aligned_reads.bam \
    | bioawk -c sam '{hist[length($seq)]++} END {for (l in hist) print l, hist[l]}' \
    | sort -n -k1

The samtools command extracts the alignements in the region given with the -L option, in sam format. The -c sam option of bioawk makes the $seq variable (corresponding to the read sequence field) available in the awk command.
The final sort ensures the size classes are ordered numerically.
We obtain a two-column size    count table.
Examples with timings:
$ time samtools view -L /Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/protein_coding.bed results_Ribo-seq_test_1/bowtie2/mapped_C_elegans/WT_1/18-30_on_C_elegans_sorted.bam | bioawk -c sam '{histo[length($seq)]++} END {for (l in histo) print l, histo[l]}' | sort -n -k1
18  201
19  560
20  1017
21  2568
22  5554
23  11221
24  27637
25  38797
26  17042
27  9853
28  15181
29  5888
30  528

real    0m0.592s
user    0m0.816s
sys     0m0.036s

Since sam format is actually simple to parse, one can try to use mawk instead for slightly higher speed (the sequence being in the 10th field):
samtools view -L regions.bed aligned_reads.bam \
    | mawk '{hist[length($10)]++} END {for (l in hist) print l"\t"hist[l]}' \
    | sort -n -k1

(Contrary to bioawk, the default output field separator is not a tabulation.)
Example:
$ time samtools view -L /Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/protein_coding.bed results_Ribo-seq_test_1/bowtie2/mapped_C_elegans/WT_1/18-30_on_C_elegans_sorted.bam | mawk '{histo[length($10)]++} END {for (l in histo) print l"\t"histo[l]}' | sort -n -k1
18  201
19  560
20  1017
21  2568
22  5554
23  11221
24  27637
25  38797
26  17042
27  9853
28  15181
29  5888
30  528

real    0m0.362s
user    0m0.500s
sys     0m0.004s

Via BEDOPS bedmap and bam2bed, along with the hist functionality in bashplotlib
$ bedmap --echo-map-size --multidelim '\n' intervals.bed <(bam2bed < reads.bam) | hist

 77|           o
 73|           o
 69|           o
 65|           o
 61|           o
 57|           o
 53|           o
 49|           o
 45|           o
 41|           o
 37|           o
 33|           o
 29|           o
 25|           o
 21|           o
 17|          oo
 13|         ooo
  9|         ooo
  5|         ooo
  1| o    o oooo
    -----------

------------------------------
|          Summary           |
------------------------------
|     observations: 118      |
|    min value: 16.000000    |
|      mean : 34.415254      |
|    max value: 36.000000    |
------------------------------

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
How can I extract normalized read count values from DESeq2 results?
The results obtained by running the results command from DESeq2 contain a "baseMean" column, which I assume is the mean across samples of the normalized counts for a given gene.
How can I access the normalized counts proper?
I tried the following (continuing with the example used here):
> dds <- DESeqDataSetFromMatrix(countData = counts_data, colData = col_data, design = ~ geno_treat)
> dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
> res <- results(dds, contrast=c("geno_treat", "prg1_HS30", "prg1_RT"))

Here is what I have for the first gene:
> res["WBGene00000001",]$baseMean
[1] 181.7862
> mean(assays(dds)$mu["WBGene00000001",])
[1] 231.4634
> mean(assays(dds)$counts["WBGene00000001",])
[1] 232.0833

assays(dds)$counts corresponds to the raw counts. assays(dds)$mu seems to be a transformation of these counts approximately preserving their mean, but this mean is very different from the "baseMean" value, so these are likely not the normalized values.

The normalized counts themselves can be accessed with counts(dds, normalized=T).
Now as to what the baseMean actually means, that will depend upon whether an "expanded model matrix" is in use or not. Given your previous question, we can see that geno_treat has a bunch of levels, which means that expanded models are not in use. In such cases, the baseMean should be the mean of the base factor in geno_treat.

mean(assays(dds)$mu["WBGene00000001",] / sizeFactors(dds)) gives something very close to res["WBGene00000001",]$baseMean (181.9139). The vignette says that assays(dds)$mu correspond to s_j * q_ij where, according to the DESeq2 paper, s_j is the size factor for library j and q_ij is, if I understand correctly, an estimation of a quantity "proportional to the concentration of cDNA fragments from the gene in the sample". Is this q_ij what you call "base factor" ?

Not really. With an "expanded model matrix" the intercept (aka, "baseMean") should be approximately the average of the normalized counts. For non-"expanded model matrices" then one level of the group factor (the "base level") becomes the intercept, so the baseMean is then just the average of its counts. The base level is the alphanumerically first level in R by default, if you have grp = factor(c("WT", "MUT")) then the MUT group is used to calculate the baseMean (assuming no expanded model matrices).

Here, I find that the mean of all normalized values (mean(counts(dds, normalized=T)["WBGene00000001",])), is actually equal to the baseMean (res["WBGene00000001",]$baseMean). This is not the case if I compute the mean only for a given level of the geno_treat factor.

I guess it using an "expanded model matrix" then. I know Michael used to code in exactly the cases when that would happen, but perhaps he's since generalized the code, which is nice :)

It depends what you mean by “normalised”. As Devon said, the normalized = TRUE argument to the count function gives you normalised counts. However, these are “only” library-size normalised (i.e. divided by the sizeFactors(dds)).
However, as the vignette explains, downstream processing generally requires more advanced normalisation, to account for the heteroscedasticity of the counts. This is often done by simply logging the counts but this has obvious drawbacks (most trivially, what do we do with 0 counts? A workaround is to add a pseudocount but that’s problematic too).
DESeq2 offers two different methods to perform a more rigorous analysis:

rlog — a regularised log, and
vst — a variance stabilising transformation.

You’d generally use either of these for downstream analysis, not count(dds, normalized = TRUE).

to make my confusion clear let's say im going for any parametric or non parametric test between any comparison for genes ,can i take the rlog or vst ? As rlog and vst i take it for clustering ,pca  ,correlation ,wgcna .So is it possible ?

@krushnachChandra It depends on the test and the comparison. DESeq for instance does not use either rlog or vst for its test of differential expression, as these would be inappropriate.

yes that i read in the paper , for example i got a cluster and i wanted to see what if the genes are having significant difference in their expression so if i have to run say a non parametric test can i use those rlog ? values or i have to use the count normalised ones ?

@krushnachChandra You cannot use count normalised values in this scenario, nor simple rlog values, since they are not normalised for transcript length (and thus longer transcripts will show higher expression on average). You could use something like size-factor normalised TPM instead. You can also log those TPMs but for a nonparametric test it would make no difference — that’s kind of the point of non-parametric tests.

So for PCA ,clustering i can use rlog i hope that part im doing correct let me know, i read your post , so can i use the deseq2 value and convert them into TPM ,but there is one part where even deseq2 also calculate size factor "dds <- estimateSizeFactors(dds)" this i guess . So can i use your TPM conversion method on deseq2 values? or just counts .

im curious about this paper https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4986529/ Figure 6 last box plot "(Right) Mean expression of the two gene sets in each population. Error bars represent SD. CPM, counts per million" where they use cpm values i  guess these are edge R output ,so here they used that value to make comparison ,so i was wondering if it can be done with rlog values ?sorry for mundane question my concept about these are very fuzzy as im not a statistician ..


text:
What kind of "gff" format does bioawk parse?
I was wondering if I could use the gff parsing capability of bioawk to facilitate the parsing of gtf files, and I looked at the following help message:
$ bioawk -c help
bed:
    1:chrom 2:start 3:end 4:name 5:score 6:strand 7:thickstart 8:thickend 9:rgb 10:blockcount 11:blocksizes 12:blockstarts
sam:
    1:qname 2:flag 3:rname 4:pos 5:mapq 6:cigar 7:rnext 8:pnext 9:tlen 10:seq 11:qual
vcf:
    1:chrom 2:pos 3:id 4:ref 5:alt 6:qual 7:filter 8:info
gff:
    1:seqname 2:source 3:feature 4:start 5:end 6:score 7:filter 8:strand 9:group 10:attribute
fastx:
    1:name 2:seq 3:qual 4:comment

I see that there are 10 fields defined for gff parsing.
However, when I look at the pages for gtf and gff in the ensembl website (gff/gtf and gff3), all have 9 fields.
I'm curious about these "filter" and "group" fields. What are they meant for? Are they extracted from some of the columns mentioned in the above pages?

I would consider the description there a bug. The filter is actually the strand, strand is the frame, group is the attribute, and attribute does nothing. These are really meant to be the 9 columns.
Edit: There's a bug report related to this.
Edit 2: I've made a pull request to clarify this and fix the aforementioned bug report.
Edit 3: I realized that I never directly answered the title of your question (mea culpa). bioawk itself will work with gff, gff3, or gtf files. It really is just treating them as tab-separated files with named columns (this is surprisingly convenient, since it's a PITA to remember what column does what).
Edit 4: The PR has been merged. If you install from github then you'll see corrected field names.


text:
Customizing bigWig file
I generate bigWig files using bamCoverage from deeptools, in part for my colleagues to visualize their mapped libraries in the IGV viewer.
A problem is that the displayed track name is apparently the file name, which is not convenient because some files for different libraries have the same name. The reason is that the files are organized with one directory per library. I would prefer to avoid changing the way I name and organize the files.
Is there a way to edit my bigWig files in order to have control of the track name independently from the file name?
As far as I understand, bigWig format derives from wiggle format, which have a track definition line that can contain a name and a definition attribute: https://genome.ucsc.edu/goldenpath/help/customTrack.html#TRACK.
Maybe there is a way to convert a bigWig file back to wiggle, edit the track definition, and then rebuild an updated bigWig file.

There's no equivalent to the wiggle header in bigWig (or bigBed) files, which is why UCSC uses the file name. This is actually the reason for the track line stuff that you linked to, since you can then specify a name and just point to where the bigWig (or other format) file is on the internet.
BTW, you can certainly convert your bigWig to wiggle, add the track line, and convert it back with UCSC tools (e.g., bigWigToWig and wigToBigWig). It shouldn't change the track name when you display it in UCSC (though if it does, then let me know, since I'll modify pyBigWig and libBigWig to support that then).

Using bigWigToWig, I can indeed obtain a wig file from the bigwig and add a custom track name that is displayed in IGV. Then, rebuilding a bigwig using wigToBigWig indeed looses the track name and I get the file name instead in IGV. This is for biwig files obtained directly with bamCoverage. I failed at reconverting wig to bigwig when the initial bigwig file was "manually" generated using pyBigWig. wigToBigWig complained as follows: line 1507430 of tests/tmp.wig: chromosome I has 15072434 bases, but item ends at 15072440. I guess I should post another question regarding this issue.

Thanks for the update, at least my understanding of bigWig headers was correct :)


text:
Dealing with absence of coverage
In my snakemake workflow, I use deeptools bamCoverage (in a make_bigwig rule), then computeMatrix and plotProfile (in different plot_*_profile rules) to create coverage profiles for various kinds of small RNAs. It sometimes happens that there are very few reads in a particular category of small RNAs. This results in bamCoverage exiting with an error (when I use the --skipNAs option, in order to make coverage profiles not take into account empty regions):
Error: The generated bedGraphFile was empty. Please adjust
your deepTools settings and check your input files.

The fact that bamCoverage fails make things a little annoying for me, because the later the failure, the easier to handle things in my workflow: I will need to catch the error and generate a dummy output file in the make_bigwig rule, and then handle the special dummy file in the plot_*_profile rules.
I would be fine having an empty coverage file: the absence of reads (either in the bam file, or matching selection criteria in bamCoverage) is not due to an error.
Is this just not possible, by the format definition?
Did I miss an option to have bamCoverage accept making empty files?
Would this make computeMatrix fail?

I'll fix the bamCoverage issue in the 2.6.0 release (it may already be fixed since I rewrite how bigWig files are made for that release). You can follow my status on the issue page. In principle you should be able to do this, it's just a deepTools limitation.
This should work fine with computeMatrix, though keep in mind that you're just going to get a big block of NA values.
Update: This error ends up coming from pyBigWig. I've just found that bedGraphToBigWig also can't make a bigWig file with now entries, so I suspect this is sort of undefined behavior. Regardless, I'll alter pyBigWig to produce a file that it can open and use in a meaningful manner.