3

What you are looking for is "Differential accessibility analysis", you are using ATAC-seq, hence "accessibility". You will need to "call" and then "normalize" "count"(*) peaks, after these, you can use DESeq2 or edgeR for the "differential" part. Here is a nice Nextflow workflow that addresses all of the steps of a differential accessibility analysis of ...


3

I have provided a solution in Biostars: If you edit the file findKnownMotifs.pl and move the call to printMotif to come before testing $svgFlag you will find that .motif files are printed again. In the following I have commented out the call from where it need not be, and introduced it as the first line, right before the if block: printMotif($_, $...


3

The answer is in the supplementary methods page 3. Finally, all fragments in the same library with duplicate start and end coordinates were removed using Picard


2

The steps you describe are correct. For step 2 it is usually normalized to mean 0 and variance 1. However the "machine learning" part is important. Having several samples being technical replicates will make the integration task easier. However, you have too few samples to make any good prediction. At most I would describe it as an exploratory analysis. I ...


2

Supplementary File 5 and 6 contain the code you're looking for.


2

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


2

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&...


2

You can see from your table that the different peaks correspond to different promoters for different transcripts, for example, FAM41C has ENST00000432963.1, ENST00000446136.1, ENST00000635557.1. And the fact that they are all within the TSS (see distanceToTSS) suggest these peaks are actually falling onto actually TSS. The underlying issue is are these TSS ...


2

You have two different peaks at different positions with different annotations. This makes complete sense, you're just misreading it.


1

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(...


1

Its highly unlikely that you will be able to categorize promoters on this basis using only ATAC data. This is a consequence of people using the terms "promoters" and "transcription start sites" interchangeably, when they are actaully different things. CAGE data measures transcription start sites which is the specific location where ...


1

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


1

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


1

I'll consider all peaks around promoter of particular gene as one peak. Moreover I'll interpret this as an open chromatin region = active gene. Why? Why there's no one wide peak which represents open chromatin in promoter of particular gene? Why is divided into several narrower peaks? There might be some transcription factors or other proteins binded to DNA, ...


1

The way we do this (and I think its quite common?), is to merge all the files of the same condition (taking the same number of reads from each), and then calling peaks on the merged sample - so you would have two peak sets - one from normal and one from disease. You then do the union of the two peak sets. This makes annotation much easier. For downstream ...


1

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


1

I would be very hesitant to use SES normalization in a case like this. We recommend that people run plotFingerprint first and see if there is good separation between the samples and only then using SES. I would not expect there to be drastic separation between the samples for ATAC-seq data, so I don't think it will work well. Our standard ATAC-seq ...


1

import numpy as np np.random.seed(0) import pyranges as pr gr = pr.random() gr.Score = np.random.randint(100, size=len(gr)) gr = gr.slack(25) # make data wider for this example print(gr) t1 = gr.tile(50) def increase_by_25(df): df = df.copy() df.Start += 25 df.End += 25 return df t2 = t1.apply(increase_by_25) tiled = pr.concat([t1, t2])...


1

The relevant help page of homer seems to be this one. On section 13 it says: An HTML page is created in the output directory named homerResults.html along with a directory named "homerResults/" that contains all of the image and other support files to create the page. These pages are explicitly created by running a subprogram called "compareMotifs.pl". ...


Only top voted, non community-wiki answers of a minimum length are eligible