jscaber
commented
2 years ago Thanks Ian!
On a related, but unaffected note: This reminds me to make a note that the 2 pass STAR method is currently only partly implemented in the mapping pipeline. The pipeline currently uses the junction file of each individual run, and I have now corrected this for the second run to use a joint junctions file from all junctions found with the first pass (with some minor filtering of low confidence junctions built in). I will make a PR very soon.
The buildJunctions function in the mapping pipeline builds a bed4 file of junctions: contig, start, end, strand
to pass to various mappers. It does this by taking a bundle of gffs from a transcript iterator, sorting them into start order, and then walking along, creating intervals from the end of one entry to the start of the next.
https://github.com/cgat-developers/cgat-flow/blob/87fcdbc71be79c8964a2af243417becc787ecea2/cgatpipelines/tools/pipeline_mapping.py#L619-L635
Unfortunately this fails to take account of the fact that the gffs list will contain things other than exons (it will include both CDS and transcript entries). In particular the transcript entries will often become the first entry when sorted on starts, but their end is after the end of all other entries. This means you get invalid entries where the start is at a higher coordinate than the end. And the interval doesn't refer to a junction.
This file is used by mapReadsWithTopHat2 and mapReadsWithHisat. One assumes no one has used TopHat2 for years . However it leads Hisat2 to produce and invalid BAM file. See:
https://github.com/DaehwanKimLab/hisat2/issues/365
PR incoming.