I have some sequencing data from a captured region that is a known paralog edited. For now, I have been mapping the data using standard minimap2 flags for PacBio DNA sequencing:

minimap2 -ax map-pb GRCh37.fa fastqfile.fastq > samfile.sam


I captured a 10kb region and what I am seeing in IGV is puzzling. The region directly before and directly after the paralogues are covered at about 1000X, but the genes themselves are not covered at all:

The two genes in the middle are the paralogues. I left them off to protect data. This data is also only the sorted data, the duplicate marked data has lower coverage (600X).

Question 1

I was wondering if there is a different approach to mapping I could take that would output all alignments for a read, and then specifically only multi-mapped reads?

Question 2

Is this process normally accomplished by minimap2, or are other aligners more effective for paralogous mapping?

Question 3

The interesting thing is, while my average read length was 6kb, I am seeing the reads here are about 4kb. I have another region with very high coverage where the reads spanned both paralogues and mapping was successful.

Question: Why am I only seeing reads that are about 4kb long when the average read length for this dataset is 6kb?

Any help on this would be great! Thanks and I hope I have been clear on my question,

Dennis

• Given that there are no spanning reads, I think that your igv screenshot is more likely to be caused by non-specific capture than by alignment. If it were me, I would double check capture probes. Dec 12 '18 at 1:11

Question 1

This command will get you all of supplementary alignments for the reads. This isn't exactly what you want though. You want all of the reads that have more than one mapping.

samtools view -f 2048 -h myfile_sorted.bam > supp_only.bam


This bash script returns a bam file that contains all alignments for reads that have multiple mappings. It works by only selecting reads that have supplementary alignments (2048 flag), or are not the primary alignment (256 flag), then using picard tools to filter the bam file. Using picard tools is the fastest way I have seen to filter a bam file given a list of reads.

#/bin/bash
BAMFILE=/path/to/my_file.bam
OUTFILE=/path/to/multiply_mapped.bam
PICARD_JAR=/path/to/picard.jar
TEMP=multiply_mapped.txt

for filt_int in 256 2048; \
do samtools view -f "$$filtint"$${BAMFILE} | \
grep -v '#' | cut -f1; done | sort | uniq > $${TEMP}; \ java -jar$${PICARD_JAR} FilterSamReads \
I=$${BAMFILE} \ O=$${OUTFILE} \
READ_LIST_FILE=$${TEMP} \ FILTER=includeReadList rm$${TEMP}


This command took 21 seconds to execute for my bam file with ~300,000 long reads mapped and several repetitive regions.

Question 2

Filtering the bam file should be pretty standard, as addressed above in my response to Q1. However, different long-read aligners are more sensitive to distantly-related paralogs. You might consider trying the blasr aligner since it is tuned for PacBio long reads, or try bwa-mem to compare your results with minimap2. Realistically though, most of these aligners will probably behave similarly to one another.

The fact that there are paralogs does not matter to the aligners. The aligner just sees these regions as two different places to which a single read can map just about equally well. When there are multiple similar regions to which a read can map, as in the case of paralogs, the aligner will report both alignments. I do not know of any long-read paralog-aware aligners.

As to why none of the reads mapped to the paralogs: Perhaps they have low mappability, as in they are repetitive? Also, maybe the capture didn't work properly? It looks like there is serious sequencing bias to the dataset given the uneven coverage.

Question 3

Are the reads that are mapping actually 4kb, or are the alignments only 4kb long? You may have many 6kb-long sequences that have been clipped at either the 5' or 3' end. If this is not a heavily-curated genome then the assembly may be incorrect around the paralogs, preventing proper mapping.

One suggestion is to take whole-genome long reads from the same organism and to map them. Then check if the mapping to this region is equally strange compared to your capture data. This will give you a better idea of what caused the uneven coverage and seeming lack of coverage around the paralogs.

• I will definitely give this a try! Thanks for the input. I will upvote this for sure but am looking for a more complete answer for handling these paralogues. I am giving blasr a try and trying some titrations of input variables. It has some advanced parameters for largely identical sequences (which is what these are). Thanks! Dec 11 '18 at 22:17
• Thanks. If you are able to edit your question to be more clear about how you would like to handle the paralogs we can try again. Right now from your question this appears to be the only sentence that says what output you would like: "I was wondering if there is a different approach to mapping I could take that would output all alignments for a read, and then specifically only multi-mapped reads?" Dec 11 '18 at 22:20
• Sorry for the ambiguity, I have updated the question to address the idea of different mappers or approaches. Thanks Dec 11 '18 at 22:23
• Great, complete answer. Thanks for the edits and update! Dec 11 '18 at 22:49

I have found an alternative solution that has been excellent for this paralog mapping issue, so I thought I'd post what I had worked on:

Question 1+2

I was wondering if there is a different approach to mapping I could take that would output all alignments for a read, and then specifically only multi-mapped reads?

Using blasr, pacbio specified some options in their github docs for long, nearly identical reads:

blasr input.bam /path/to/reference.fa --hitpolicy all --bam --out alignments.bam --minMatch 8 --maxMatch 15 --nproc 4


The link to the docs as to why this is useful is here, although I did tweak the parameters to smaller seed sizes to create more anchors to the reference, which produced more secondary alignments. This is effective because this will output up to 10 secondary alignments per read, which I saw quite frequently for my region of interest.

After this, sorting was accomplished using Picard SortSam, with the following flags:

java -jar picard.jar SortSam INPUT=$${SAMPLE} OUTPUT=./sorted/$${SAMPLE}\
SORT_ORDER=coordinate VALIDATION_STRINGENCY=LENIENT\
CREATE_INDEX=true MAX_RECORDS_IN_RAM=100000


The blasr files have much more information than normal sams produced by aligners like minimap2 and bwa mem, so the MAX_RECORDS_IN_RAM=100000 as opposed to 500000 helps to ensure less frequent job failures (this was done on a cluster).

Next, I parsed out the region of interest:

samtools view -h input.bam "chrX:1000-2000" > extracted.bam # this is fake


awk '$5==0' extracted.bam > mapq0.txt  Get read names to grep for later: awk '{print$1}' mapq0.txt > read_names_mapq0.txt


Search for all reads with readname from the original bam to get mapping positions:

samtools view infile.bam | grep -f read_names_mapq0.txt > all_reads_mapq0_mapq_high.txt


Keep only columns: 1) read_name, 2) Sam Flag, 3) Chr, 4) Pos, 5) MapQ:

awk '{print $$1,$$2,$$3,$$4,\$5}' OFS='\t' all_reads_mapq0_mapq_high.txt > mapq_read_info.txt


Then I wrote a simple python script to sort by read name, count the number of times a read shows up, and see all the positions it maps to. An example of the output:

        Read_Name       Flag    Chr     Pos     MAPQ
...


        Read_Name       Flag    Chr     Pos     MAPQ

Flag=0: read was mapped to + strand

• 100% identity doesn't count if the alignment only covers a small fraction of the query. Those fragmented alignments are probably hits to repeats. Minimap2 can be tuned to output such hits with -p.1 or -P, but most of extra alignments are not useful. Dec 15 '18 at 21:30