10

The key function is call_broadpeaks: The description attached to the function says: This function try to find enriched regions within which, scores are continuously higher than a given cutoff for level 1, and link them using the gap above level 2 cutoff with a maximum length of lvl2_max_gap. scoring_function_s: symbols of functions to ...


9

Aside: Cross-correlation is largely meaningless, regardless of what some of the ENCODE folks might argue. When we process our DEEP samples we don't even look at that value. Regardless, if you're using SPP/phantomPeakQual for cross-correlation then note that it already removes the highest peaks from your dataset before computing the cross-correlation (in ...


7

I used this in the past for ChIP-seq data and it generated SNVs: samtools mpileup \ --uncompressed --max-depth 10000 --min-MQ 20 --ignore-RG --skip-indels \ --fasta-ref ref.fa file.bam \ | bcftools call --consensus-caller \ > out.vcf This was samtools 1.3 in case that makes a difference.


7

I don’t think it matters. Both are easy to merge (BAM via samtools merge, and (gzipped) FASTQ via cat), and neither method has specific disadvantages, unless your FASTQ files are sorted for some reason (but they generally shouldn’t be). One advantage of keeping the FASTQ files separate is that it makes it slightly easier to parallelise the mapping step: ...


5

Another approach is htsbox. You can get a candidate list with: htsbox pileup -Cvcf ref.fa -q20 -Q20 -s5 file.bam > out.vcf Here, -q sets min mapping quality, -Q sets min base quality, -v outputs variants only -c outputs VCF, -C gives you base counts on both strands and finally -s5 requires at least 5 high-quality bases to call out an allele. It is ...


5

MSigDB has a collection (C3:TFT) of gene sets corresponding to transcription factor targets. Harmonizome has functional terms for genes extracted from over a hundred publicly available resources.


4

Your original command without --nomodel --extsize ... is probably the most accurate. This warning stems from a time when reads were much much shorter and likely never made that much sense to begin with. Broad peak calling in MACS2 basically works by finding a bunch of nearish narrow peaks and merging them. If you have really broad signals then use something ...


4

I strongly suggest that you not try to come up with your own package for this when things like CSAW already exist in bioconductor and provide a number of useful normalization options. For visualization purposes, I think it's best to simply take the log2 ratio (e.g., using bamCompare from deepTools) or alternatively normalize by total coverage (possibly ...


4

That's a fastq file, you will want to align it to the genome, call peaks, and then use something like MEME to determine binding motifs.


4

In general, there is two options to identify targets for transcription factors: experimental (ChIP-seq) and sequence-based predictions. TF binding from experimental data The are multiple projects that produce binding data of transcription factors and quantify their peaks across the genome. The advantage here is that you know that binding actually occurs, ...


4

iRegulon takes a sequence-based approach to finding transcription factor targets. There's a Cytoscape app that you can use to find the regulators of a given gene list, or the targets of a particular transcription factor. Transcription factor binding sites are predicted using a collection of position weight matrices (PWMs) from a number of sources, including ...


4

I would ignore peak calling for this and instead compute enrichment of ChIP/input for the genome (e.g., with deepTools or presumably homer) and then plot it for the genes of interest individually (e.g., using IGV or pyGenomeTracks) or as a group (e.g., with computeMatrix). If the peaks are obvious and you trust your peak calling then sure you can just use ...


4

There is no such thing as a hypergeometric test, at least in statistical textbooks. It's a fisher test based on hypergeometric distribution. If it is chip-seq for the same target, i.e biological replicates, significance of overlap is not quite meaningful. You get more information by for example checking the correlation of the coverage between your ...


3

This paper https://www.biorxiv.org/content/10.1101/728519v1.full examines the mapping of some human telomeric sequences, and indicates that there is something like 130kb missing at the telomeres in the HG38 assembly. One option might be to take the HG38 assembly, and add some missing telomere sequences as an additional pseudochromosome to the FASTA file. ...


3

You could make your own such file with FIMO, the JASPAR MEME files available via the MEME download site, and your genome build of interest (per-chromosome FASTA files from UCSC goldenpath, say), e.g., for assembly hg19 and UCSC chromosome naming scheme: $ for chr in `seq 1 22` X Y; do echo ${chr}; wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/hg19/...


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

I'd say mostly this is a question of understanding the underlying biology and the relevant literature. If it is not known in the literature whether a mark is peaky or broad, evidence might come from FISH studies or Low throughput qPCR. Another way to look would be to examine the signal expressed as fold enrichment over input on a genome browser and look if ...


3

For ChIP-seq it shouldn't really matter. But do be aware that by default, samtools merge retains read group information (the @RG field in the header) from each input file. This could pose a problem for some downstream analyses (e.g. for the GATK HaplotypeCaller) if you want the merged data to be considered as all part of the same sample. You can change this ...


2

Agree with the others that it doesn't really matter. One thing to note though - if you're deduplicating your BAM files (you probably should for ChIP-seq data), make sure that you do this after merging.. :)


2

This is very good question, but as far as I know there really isn't a good answer. I will attempt to provide some comments and tips. However, I'm trying to do similar with a dataset which may consist of repeats (I have previously aligned the dataset using bowtie, and a lot of reads aligned to repeats). If this dataset is also PE and the alignment is ...


2

I think you might have changed the separator (or at least have some kind of inconsistency from the required format) for your file. Note that peak output files from MACS2 are variants of BED files. It seems you need to have a tab separator for this type of file. I copied your example file and then ran awk -v OFS="\t" '$1=$1' your_peaks.narrowPeaks > ...


2

You can use Homer to retrieve tag densities, quoting from their manual: Calculating ChIP-Seq Tag Densities across different experiments annotatePeaks.pl is [a] useful program for cross-referencing data from multiple experiments. In order to count the number of tags from different sequencing experiments, you must first create tag directories ...


2

You could try setting a high p-value threshold when you're calling peaks to retain the "non-significant" values. Something like: macs2 callpeak -t <ChIP>.bam -c <Control>.bam -f BED -g hs -n test -B -p 0.7 You probably don't want to set the threshold to 1, since you'll likely want to avoid some noise. But you can then use these "non-...


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

This question is somewhat generic, so a generic answer is that ENCODE has a Transcription Factor ChIP-seq Data Standards and Processing page that can give you a useful starting point. For TF ChIP-seq data with replicates, the Irreproducible Discovery Rate (IDR) method helps leverage replicates to produce higher confidence peak calls, producing both "optimal"...


2

You could first look at the degree of correlation between the two replicates - what proportion of peaks are shared between the two, versus peaks found in one sample only? This will give an idea of how repeatable the analysis is, and how many peaks are variable between samples or due to artefacts of the method. When it comes to interpreting the biology ...


2

bedtools intersect with multiple files after -b performs pairwise intersections between -a and each file in -b. It does not intersect all files simultaneously. See here, from the Bedtools documentation. The green section only contains intervals from -a that also intersect at least one interval in any -b file. What I think you're looking for is bedtools ...


2

Headers should not contain tabs: track type=narrowPeak name="Somite narrowPeak" description="Somite narrowPeak" Ensure that you have NO tabs on that line.


2

Depends how the files were made. In the simplest case yes, the height represents the pipeup of reads from the BAM file that was used. For a direct (visual) comparison you have to normalize the files though as otherwise sequencing depth confounds the height of the peaks, like if file A has ten times more reads and the peak is 10 times higher than in B then ...


1

Before you proceed, it is usual to do quality control on the .fastq data. The FASTQC program is a good place to start. Once you have checked the quality, trimming is the typical next step to remove low-quality data. There are options such as Trimomatic, trimGalore for this. Once you are sure you have clean data, your first step is alignment - align the ...


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