2
MRNM stands for "Mate reference index". So Picard found something in the RNEXT field which should be set only for paired-end reads but the rest of the file looks like single-end.
The problematic line in your code is:
line.next_reference_id = 0
This sets the RNEXT SAM field to whatever Pysam stores as a reference with index 0 (next_reference_id). ...
2
Like Devon pointed out, most likely you should sort out whether the files have been marked for duplicates correctly.
You can also use samtools rmdup too
samtools rmdup <bamfile> <output>
Also, I find it very odd that 50% of the reads are lost when you keep only properly paired reads.
I would suggest removing the duplicates and then ...
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
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
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
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
In my opinion, the easiest way to merge bed files is to use bedtools merge.
1
Convert each of your Excel spreadsheets to tab-delimited text files (how-to).
Install and run Cygwin, if using Windows. If you are using OS X, open up the Terminal app.
Install the BEDOPS toolkit to get sorting and set operation tools (if you're not averse to using them).
Strip out the Microsoft line endings, strip the header, and sort each of the N files:
$...
1
You can use rtracklayer from R bioconductor to import it, filter and do much more stuff, i use the code from this blog for importing the narrow bed file. Below I use an example file since you did not provide yours:
library(rtracklayer)
fl = "zas1_default_peaks.narrowPeak"
extraCols_narrowPeak <- c(signalValue = "numeric", pValue = &...
1
You can do it by awk.
awk -v OFS="\t" '{ if ( ($3-$2) >= 500 && ($3-$2) <= 1000) { print } }' in.bed > out.bed
Assuming there is no header.
1
For sharing, you can test the proportion of overlapping peaks using bedtools intersect or the find.overlap function in GenomicRanges packages in R.
What makes more sense is perhaps to merge all the peaks together using bedtools merge, and then count the number of reads/fragments in each peak for each sample using featureCounts. The output of this is a matrix,...
1
Venn diagrams can be helpful, but for peak calling they can be somewhat misleading.
This topic is addressed in the DiffBind documentation, in the "Comparison of occupancy and affinity based analyses" section.
Suffice to say, asking "which peaks are shared?" is not the same as asking "is there evidence for differential peaks, and where are they?".
For this ...
1
Check out the tool Intervene. It's exactly designed for this. Github link.
1
not gnuparallel but nextflow. I don't know macs and what is 'first.bam' and 'second.bam' is not clear not me but I hope you'll get the idea.
The nextflow:
Channel.fromPath(params.bams).splitCsv(header: false,sep:'\t',strip:true).map{T->file(T[0])}.into{bams1;bams2}
boolean isControl(file) {
return file.name.endsWith("_Input_sorted.bam_rem.bam") || ...
1
You need to provide sorted bam files, ideally with an index. If you run ChIPQCsample, it calls ChIPQC:::sampleQC which iterates through the chromosomes or contigs, using the index, and the absence of the index will crash it.
We can use an example, first in linux:
wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibTfbs/...
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])...
Only top voted, non community-wiki answers of a minimum length are eligible