Hello Pierre,

Your script is very useful in my research. And is there a way to dissect with multiple-base pattern in a specified region from BAM file? (Like extract alignment with chr19:100000-100003 as ATCG in reverse strand). Another question is I noticed that in the manual of samjdk there is a specification: --pair [20171110] PAIR-MODE .The signature of java function is public Object apply(final List<SAMRecord> records). This function must return true to accept the whole list, false to reject eveything, or another List<SAMRecord>.Input MUST be sorted on query name using picard SortSam (not samtools sort https://github.com/samtools/hts-specs/issues/5 ). Does that mean when I using pairend BAM files as input, I should use picard to sort them first?

Thank you!

Hi Pierre,

Thank you for the script. It has helped me a lot to select reads that contain a substitution mutation in my cDNA sequencing data. However, I would like to do the same thing for selecting reads that have indel mutations (both 1 bp and multiple bp), and splice site mutations. The splice site mutations are located intronically, and thus not always present in my reads, so a selection for these would have to be based on splicing (such as exon skipping). Do you know a way to do this using your program?

Thank you!

Hi Pierre, Thank you very much for making this tool. Every time I need a not-so-standard tool there is one piece of software in jvarkit that does specifically what I want. This time I am interested in selecting reads supporting one allele, but I'd like to know how to select several positions specifying one base each time. I've tried this script.js, based on a simple modification of your script.js:

final String contig = "scaffold_3";
final int[] mutpos = {9379104, 9379106, 9379131, 9384562, 9414473, 9416683, 9416773, 9417122, 9417182, 9417260};
final char[] mutbase = {'T', 'G', 'A', 'T', 'T', 'A', 'A', 'G', 'C', 'C' };
int arrayLength = mutbase.length;
if(record.getReadUnmappedFlag()) return false;
if(!record.getContig().equals(contig)) return false;
for (int i = 0; i < arrayLength; i++){
  if(record.getEnd() < mutpos[i]) return false;
  if(record.getStart() > mutpos[i]) return false;
  int readpos = record.getReadPositionAtReferencePosition(mutpos[i]);
  if(readpos<1) return false;
  readpos--;
  final byte[]    bases= record.getReadBases();
  if(bases[readpos]==mutbase[i]) return true;
  return false;
}
return false;

However, this ends up giving me less entries than using a single position.

Please note that this is my very first attempt to "code" anything in java. This does run, which I did not expect, but the results are not what I want. Is there a way to ask samjdk to return a read whenever a mutpos matches its corresponding mutbase?. I don't care if I got duplicated reads, I can sort and uniq them later (I'm testing some hundreds of SNPs, so the sort should be rather fast).

Of course a different option would be to run this once per SNP in a loop (that is, generate a script.js per SNP), but I'd rather run all at once.

Thank you again for writing jvarkit.

EDIT: Apparently the issue came from the return statement, that ends the for loop and finishes the script. I've modified the script and it works, I'm just leaving it here for anyone who comes in the future (maybe even myself!)

final String contig = "Your_Contig";
final int[] mutpos = {Pos1, Pos2, Pos3, PosN};
final char[] mutbase = {'BaseForPos1', 'BaseForPos2', 'BaseForPos3', 'BaseForPosN' };
int arrayLength = mutbase.length;
if(record.getReadUnmappedFlag())  return false;
if(!record.getContig().equals(contig)) return false;
ArrayList<htsjdk.samtools.SAMRecord> ret = new ArrayList<htsjdk.samtools.SAMRecord>();
for (int i=0; i < arrayLength; i++){
  if(record.getEnd() < mutpos[i]) continue;
  if(record.getStart() > mutpos[i]) continue;
  int readpos = record.getReadPositionAtReferencePosition(mutpos[i]);
  if(readpos<1) continue;
  readpos--;
  final byte[]    bases= record.getReadBases();
  if(bases[readpos]==mutbase[i]) {
    ret.add(record);
  }
  continue;
}
return ret;

I hope this is useful for someone else.