domainDirections

[francisfa]

Hi,
I want to compare the different boundary strength between 2 samples. So I using the Chapter 8 code like:
finder <- domainDirections(input, hm.param, width=1e5, span=10) to calculate the DI first.

But, there is the error:
Error in preloader(files, param=param, retain=retainer): 'intervals' in 'HiCXXX.h5' is not present in fragment chromosomes

[LTLA]

Best guess is that your file has a chromosome that isn't present in hm.param$fragments. In particular, intervals doesn't sound like it's a real chromosome name; you should probably check the originating BAM file and see if it's present there.

[VonAlphaBisZulu]

Unfortunately, even when you use the same genome for alignment and fractioning, diffHiC may still throw this issue. This happens when a certain chromosome does not have any restriction site, but you have alignment to this chromosomes in your bam.

Your only option is then to remove all alignments to these chromosomes from your bam and from the bam header and try again. I personally think diffHiC should work in that case, regardless. It deters people from using the package.

Here is some example code. Let's assume you have (uniquely and primary) aligned paired-end reads, generated separately with bowtie2, bwa mem or similar and stored the alignments in aligned1.bam and aligned2.bam (again, no multimappers):

# Add the forward and backward flag to the aligned reads

samtools view -bu --add-flags 65 aligned1.bam | samtools sort -u -n aligned1_sorted.bam

samtools view -bu --add-flags 129 aligned2.bam | samtools sort -u -n aligned2_sorted.bam

 

# Merge name-sorted bam files of aligned reads (alternating lines from each file)

samtools merge -f -n merged.bam aligned1_sorted.bam aligned2_sorted.bam

 

# Remove all reads whose chromosome name is longer than 7 letters or is chrMT, then remove second reads that lost their pair earlier

# Here we process the header separately from the reads. In both cases we use awk to strip the relevant info

( samtools view merged.bam -H | awk '!($1=="@SQ" && length($2)>10) && $2!="SN:chrMT"'& \

samtools view merged.bam | awk 'length($3)<7 && $3!="chrMT"' | awk '{if ($1 in used ) {if (used[$1] != "" ) print used[$1] ; print ;} else used[$1]=$0 ;}' ) | \

samtools sort -u -n filtered.bam

 

# fix mate-pair flags

samtools fixmate -p -@ filtered.bam final.bam