Hello, I would greatly appreciate some input into this issue. I put all command lines I used at the bottom.
I have an RNA IP and input dataset aligned with STAR and when I input the generated bam files into MACS2 to generate peaks, MACS2 is estimating the fragment size around 14151.2bp (actual size is around 150bp) and it gives the error: AssertionError: 1000 can't be smaller than 14151.1923828125! because the default slocal parameter is 1000. So it's taking the genomic distance into consideration and not the actual mRNA sequences (inserting the introns which gives a large fragment size), even though I set --nomodel. I aligned the same dataset with HISAT2 and I didn't have this issue, it's only occurring with STAR. Please see below the command lines I used:
for mapping with STAR:
STAR --genomeDir ~/Ref_Genome/GenomeDir_STAR --runThreadN 60 --readFilesIn IP1_R1.fq.gz IP1_R2.fq.gz --outFileNamePrefix /STAR_Align --outSAMtype BAM SortedByCoordinate --outSAMunmapped Within --readFilesCommand gunzip -c --outFilterScoreMinOverLread 0 --outFilterMatchNminOverLread 0 --outFilterMatchNmin 0
MACS2 peak calling:
macs2 callpeak -t IP.bam -c input.bam -n CNT1 --nomodel --extsize 75 -g 382e6 -f BAMPE --verbose 3 --outdir ~/Documents/Bruno/MACS2
Also tried adding: --nolambda, I tried changing the -g perameter to(hs, 100e6, 2.7e9.
I am not sure if it's a problem with the STAR alignment (which could be from the index generated or from the command line I used), or if it's from MACS2 (missing some perameters?). Others have used macs2 with STAR-aligned RNA samples very similar to mine so it should work in theory.
I appreciate any feedback. Thanks!
I have met a similar error, and the result of hisat2 shared the same error too. I am thinking about using the bwa-mem to align the RIP-seq data.