paired-end pileup with annotation file - question

[joegeorgeson]

Dear JACUSA2 team,

When using pileup with paired-end alignments and passing an annotation file, if R1 maps within one annotation but R2 maps outside (either in a different annotation or no annotation), how are R1/R2 counted? Is there a way this can be controlled?

[piechottam]

Would do you mean by annotation file? A fasta reference via "-R REF-FASTA"?
The fasta reference should be the same that was used for mapping and it influences that "ref_base" column in the JACUSA2 output. JACUSA2 will ignore all non primary alignments.

Or you mean R1 maps to chr1 and R2 maps to chr2 or is not mapped. In that case JACUAS2 pileup "will count what is there".
But you can change that behaviour by filtering with "samtools flags".

The relevant part of the source code is in lib.cli.parameter.ConditionParameter.java:

public boolean isValid(SAMRecord samRecord) {
  if (! samRecord.getReadUnmappedFlag()
		  && ! samRecord.getNotPrimaryAlignmentFlag() // ignore non-primary alignments CHECK
		  && (mapq < 0 || mapq >= getMinMAPQ()) // filter by mapping quality
		  && (getFilterFlags() == 0 || (getFilterFlags() > 0 && ((samRecord.getFlags() & getFilterFlags()) == 0)))
		  && (getRetainFlags() == 0 || (getRetainFlags() > 0 && ((samRecord.getFlags() & getRetainFlags()) > 0)))
		  && errors == null // isValid is expensive
		  ) { // only store valid records that contain mapped reads
	  // custom filter 
	  for (final MaxValueSamTagFilter samTagFilter : getSamTagFilters()) {
		  if (samTagFilter.filter(samRecord)) {
			  return false;
		  }
	  }
  
	  // no errors found
	  return true;
  }
  
  // print error messages
  if (errors != null) {
	  for (SAMValidationError error : errors) {
			  AbstractTool.getLogger().addError(error.toString());
	  }
  }
  
  // something went wrong
  return false;

[joegeorgeson]

Would do you mean by annotation file? A fasta reference via "-R REF-FASTA"?

By annotation I mean bed file, the -b <BED> parameter. A simple example describing a situation of importance for me; if geneA = chr1:1-100 & geneB = chr1:110-200 (no overlap between genes), with a paired-end read where R1 maps to chr1:50-100 (only geneA) and R2 maps from chr1:110-160 (only geneB) will these reads be counted?

[piechottam]

-b <BED> defines the search region(s) of JACUSA2.
All position (including partially overlapping reads) that overlap region(s) in will be counted.

Are you missing some reads when you use -b <BED>?

[joegeorgeson]

Are you missing some reads when you use -b <BED>?

Not that I'm aware of, I just wanted to make sure I understood how paired-end reads are counted when mapped across defined regions.

Thanks!