<p>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:</p>

<p><strong>Question 1+2</strong></p>

<p>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?</p>

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

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

<p>The link to the docs as to why this is useful is <a href="https://github.com/PacificBiosciences/blasr/wiki/To-align-long-and-almost-identical-sequences" rel="nofollow noreferrer">here</a>, 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.</p>

<p>After this, sorting was accomplished using Picard SortSam, with the following flags:</p>

<pre><code>java -jar picard.jar SortSam INPUT=<span class="math-container">${SAMPLE} OUTPUT=./sorted/$</span>{SAMPLE}\
SORT_ORDER=coordinate VALIDATION_STRINGENCY=LENIENT\
CREATE_INDEX=true MAX_RECORDS_IN_RAM=100000
</code></pre>

<p>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).</p>

<p>Next, I parsed out the region of interest:</p>

<pre><code>samtools view -h input.bam "chrX:1000-2000" &gt; extracted.bam # this is fake
</code></pre>

<p>Parse out mapq=0 reads:</p>

<pre><code>awk '$5==0' extracted.bam &gt; mapq0.txt
</code></pre>

<p>Get read names to grep for later:</p>

<pre><code>awk '{print $1}' mapq0.txt &gt; read_names_mapq0.txt
</code></pre>

<p>Search for all reads with readname from the original bam to get mapping positions:</p>

<pre><code>samtools view infile.bam | grep -f read_names_mapq0.txt &gt; all_reads_mapq0_mapq_high.txt
</code></pre>

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

<pre><code>awk '{print <span class="math-container">$1,$</span>2,<span class="math-container">$3,$</span>4,$5}' OFS='\t' all_reads_mapq0_mapq_high.txt &gt; mapq_read_info.txt 
</code></pre>

<p>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:</p>

<pre><code>        Read_Name       Flag    Chr     Pos     MAPQ
read1/10158955/ccs       272     19      3395674 0
read1/10158955/ccs       0       16      217682  254
read1/10158955/ccs       272     12      64990527        0
read1/10158955/ccs       256     16      222913  0
read1/10158955/ccs       272     19      21328060        0
read1/10158955/ccs       256     19      3224052 0
read1/10158955/ccs       272     17      34878621        0
read1/10158955/ccs       272     16      84204312        0
read1/10158955/ccs       256     1       244788680       0
read1/10158955/ccs       272     2       222951342       0
read2/10223868/ccs       16      16      217057  254
read2/10223868/ccs       256     20      3078724 0
read2/10223868/ccs       272     6       133521134       0
read2/10223868/ccs       256     7       152515722       0
read2/10223868/ccs       256     21      45258117        0
read2/10223868/ccs       256     1       249160685       0
read2/10223868/ccs       272     16      223631  0
read3/10551872/ccs       272     GL000220.1      144982  0
read3/10551872/ccs       256     18      72831954        0
...
</code></pre>

<p>read_counts output:</p>

<pre><code>        Read_Name       Flag    Chr     Pos     MAPQ
read1/10158955/ccs       10      10      10      10
read2/10223868/ccs       7       7       7       7
read3/10551872/ccs       9       9       9       9
</code></pre>

<p>The numbers in each column correspond to the number of supplementary reads seen. 
A brief explanation of sam flags which I feel is relevant for unpaired reads:</p>

<pre><code>Flag=0: read was mapped to + strand
Flag=16: read was mapped to - strand
Flag=256: read was mapped to + as secondary alignment
Flag=272: read was mapped to - as secondary alignment
</code></pre>

CC BY-SA 4.0
2018-12-15T01:33:35.870Z
Minimap2 gives you up to 5 secondary alignments by default. IGV just doesn't show them.
user172818
CC BY-SA 4.0
2018-12-15T16:43:02.210Z
Yeah I know, but I tried many varying -H and -k parameters and never saw any secondary reads to my region. I know they should be there because when I chunk my reads and blast them I see 100% identity to this region as well as a few others but minimap2 is not outputting any secondary alignments there.
d_kennetz
CC BY-SA 4.0
2018-12-15T16:44:03.457Z
And I don’t just mean in the viewer, I mean in the sam file. Unless there is more I can do?
d_kennetz
CC BY-SA 4.0
2018-12-15T21:30:21.873Z
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.
user172818
CC BY-SA 4.0
2018-12-14T23:34:41.653Z
2018-12-14T23:34:41.653Z
d_kennetz
<h3>Question 1</h3>

<p>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.</p>

<pre><code>samtools view -f 2048 -h myfile_sorted.bam &gt; supp_only.bam
</code></pre>

<p>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 <code>bam</code> file.  <a href="https://bioinformatics.stackexchange.com/a/3390/2085">Using picard tools is the fastest way I have seen to filter a <code>bam</code> file given a list of reads.</a></p>

<pre class="lang-bsh prettyprint-override"><code>#/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 "<span class="math-container">$filtint" $</span>{BAMFILE} | \
  grep -v '#' | cut -f1; done | sort | uniq &gt; <span class="math-container">${TEMP}; \
java -jar $</span>{PICARD_JAR} FilterSamReads \
  I=<span class="math-container">${BAMFILE} \
  O=$</span>{OUTFILE} \
  READ_LIST_FILE=<span class="math-container">${TEMP} \
  FILTER=includeReadList
rm $</span>{TEMP}
</code></pre>

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

<h3>Question 2</h3>

<p>Filtering the <code>bam</code> file should be pretty standard, as addressed above in my response to <strong>Q1</strong>. However, different long-read aligners are more sensitive to distantly-related paralogs. You might consider trying the <a href="https://github.com/PacificBiosciences/blasr" rel="nofollow noreferrer">blasr aligner</a> 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.</p>

<p><strong>The fact that there are paralogs does not matter to the aligners.</strong> 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. <strong>I do not know of any long-read paralog-aware aligners.</strong></p>

<p>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.</p>

<h3>Question 3</h3>

<p>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. </p>

<p>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.</p>

CC BY-SA 4.0
2018-12-11T22:17:46.853Z
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!
d_kennetz
CC BY-SA 4.0
2018-12-11T22:20:09.623Z
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?"
conchoecia
CC BY-SA 4.0
2018-12-11T22:23:10.013Z
Sorry for the ambiguity, I have updated the question to address the idea of different mappers or approaches. Thanks
d_kennetz
CC BY-SA 4.0
2018-12-11T22:49:40.800Z
Great, complete answer. Thanks for the edits and update!
d_kennetz
CC BY-SA 4.0
2018-12-11T22:08:14.657Z
2018-12-11T22:47:12.230Z
2018-12-11T22:47:12.230Z
conchoecia
<p>I have some sequencing data from a captured region that is a known paralog <strong>edited</strong>. For now, I have been mapping the data using standard minimap2 flags for PacBio DNA sequencing:</p>

<pre><code>minimap2 -ax map-pb GRCh37.fa fastqfile.fastq &gt; samfile.sam
</code></pre>

<p>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:</p>

<p><a href="https://i.sstatic.net/ecXKR.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/ecXKR.png" alt="enter image description here"></a></p>

<p>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). </p>

<h3>Question 1</h3>

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

<h3>Question 2</h3>

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

<h3>Question 3</h3>

<p>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.</p>

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

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

<p>Dennis</p>

CC BY-SA 4.0
2018-12-12T01:11:06.167Z
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.
user172818
<long-reads><pacbio><mappability><minimap2><paralogs>
Mapping Reads to Known Gene Paralogs with Long Read Technology
# Mapping Reads to Known Gene Paralogs with Long Read Technology

(Asked by: username_0 on 2018-12-11)

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

**Tags:** long-reads, pacbio, mappability, minimap2, paralogs
> **Comments:**
> **username_1**: 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.  

---

## 2 answers

### Answer by username_2 on 2018-12-11. Score: 3 | Accepted answer

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.

> **Comments:**
> **username_0**: 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!  
> **username_2**: 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?"  
> **username_0**: Sorry for the ambiguity, I have updated the question to address the idea of different mappers or approaches. Thanks  
> **username_0**: Great, complete answer. Thanks for the edits and update!  

---

### Answer by username_0 on 2018-12-14. Score: 0

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

Parse out mapq=0 reads:
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
read1/10158955/ccs       272     19      3395674 0
read1/10158955/ccs       0       16      217682  254
read1/10158955/ccs       272     12      64990527        0
read1/10158955/ccs       256     16      222913  0
read1/10158955/ccs       272     19      21328060        0
read1/10158955/ccs       256     19      3224052 0
read1/10158955/ccs       272     17      34878621        0
read1/10158955/ccs       272     16      84204312        0
read1/10158955/ccs       256     1       244788680       0
read1/10158955/ccs       272     2       222951342       0
read2/10223868/ccs       16      16      217057  254
read2/10223868/ccs       256     20      3078724 0
read2/10223868/ccs       272     6       133521134       0
read2/10223868/ccs       256     7       152515722       0
read2/10223868/ccs       256     21      45258117        0
read2/10223868/ccs       256     1       249160685       0
read2/10223868/ccs       272     16      223631  0
read3/10551872/ccs       272     GL000220.1      144982  0
read3/10551872/ccs       256     18      72831954        0
...

read_counts output:
        Read_Name       Flag    Chr     Pos     MAPQ
read1/10158955/ccs       10      10      10      10
read2/10223868/ccs       7       7       7       7
read3/10551872/ccs       9       9       9       9

The numbers in each column correspond to the number of supplementary reads seen. 
A brief explanation of sam flags which I feel is relevant for unpaired reads:
Flag=0: read was mapped to + strand
Flag=16: read was mapped to - strand
Flag=256: read was mapped to + as secondary alignment
Flag=272: read was mapped to - as secondary alignment

> **Comments:**
> **username_1**: Minimap2 gives you up to 5 secondary alignments by default. IGV just doesn't show them.  
> **username_0**: Yeah I know, but I tried many varying -H and -k parameters and never saw any secondary reads to my region. I know they should be there because when I chunk my reads and blast them I see 100% identity to this region as well as a few others but minimap2 is not outputting any secondary alignments there.  
> **username_0**: And I don’t just mean in the viewer, I mean in the sam file. Unless there is more I can do?  
> **username_1**: 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.  

---

