Align paralogous contigs to reference

If you applied the --keep-paralogs flag in the SECAPR find_target_contigs function, the function will print a text file with paralogous information into the subfolder of each sample. This file contains the information about the locus id (first column) in the reference file, followed by a list of contig headers that were matched to this locus:



In [8]:

    
%%bash
head -n 10 ../../data/processed/target_contigs_paralogs/1061/info_paralogous_loci.txt









    



214	86725	263626																												
213	107618	194246																												
669	53180	125308																												
762	122510	147847	263530																											
763	153094	19254	19255	114710	177616	211218	256601																							
761	5580	46418	249536																											
663	114228	114230																												
662	33701	162610																												
764	18002	46660	184621																											
765	48851	171788	217113	218255	254469

We can use the SECAPR paralogs_to_ref function to extract the sequences of all these contigs and align them with their respective reference sequence. This will give an idea of where the contigs map on the reference and can help decide if these are truly potential paralogous or more likely the result of non-optimal contig assembly (clusters of homologous reads are broken into separate contigs) or if contigs from other parts of the genome map to the reference due to the presence of e.g. repetitive regions or other common sequence patterns.

We need to provide the function the following input items:

path to the de novo contig files
path to the reference fasta file, which was used to extract target contigs
path to the extracted target contigs (output of find_target_contigs function)

The command looks as follows:

secapr paralogs_to_ref --contigs ../../data/processed/contigs --reference ../../data/raw/palm_reference_sequences.fasta --target_contigs ../../data/processed/target_contigs_paralogs_info --output ../../data/processed/paralogs_to_reference

Depending on how many paralogous loci were identified in your samples this can take several minutes. The script will store the final alignments of the contigs and the reference sequences for each sample in the paralog_alignments folder in the provided output path. Let's look at one exemplarly alignment. You can view alignments using alignment viewers such as e.g. AliView:



In [10]:

    
from IPython.display import Image, display
img1 = Image("../../images/paralog_contig_alignment.png",width=1000)
display(img1)

In this example you see the reference sequence in the first line starts with a repetative sequence (CTCTCTCTC....). Many contigs that only consist of this repetitive sequence appear to map to it. Only one of the contigs, however, apears to be a true match, as it exceeds the repetative region and has a lot of overlap with the reference sequence. A locus like this doesn't need to be excluded, since we see that even though it's flagged as paralogous, there are no true signs of paralogy.

Back