# Genome Alignment with LAST
By Tanya Lama

**Objective**:
Write a functional script for eventually mapping our 58 low-coverage Canada lynx genomes to either the annotated Canada lynx reference genome or some other well-assembled and thoroughly annotated reference genome for downstream analysis. Here we practice alignment, explore a couple of the available Felid genomes available on NCBI, and examine output files and alignment quality scores.

**Requirements**:
For aligning to a mammal genome, you'll need a few dozen gigabytes of memory.
First, get a reference genome sequence, in FASTA format and install the latest version of LAST (>= 979).

### Part 1: Prepare a genome with pre-alignment repeat-masking
#We're masking repeat regions here as a means of increasing mapping speed and reducing compute time and memory. Features that map to repeats are often uninformative, so masking repeats is appropriate in some instances and can make things orders-of-magnitude faster without losing information. However, repeat-masking also harms alignment accuracy (by hiding some correct alignments). Here we have DNA reads with multiple coverage of a mammal genome, so masking repeats is probably a okay. We'll be more thoughtful about this and will likely choose a slower and more accurate approach when aligning our low-coverage whole genomes (this is achievable with different parameterization in *LAST*)

*LAST* can detect simple repeats such as atatatatatatat, but not interspersed repeats: for that we use *windowmasker*.

#Download and install the standalone [NCBI Blast setup for UNIX](ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/
) which includes *windowmasker*

cd /Users/tanyalama/ncbi-blast-2.8.1+-src/c++
./configure
cd /Users/tanyalama/ncbi-blast-2.8.1+-src/c++/ReleaseMT/build
make all_r

#The compiled executables will be found in /Users/tanyalama/ncbi-blast-2.8.1+-src/c++/ReleaseMT/bin

#Consult the *windowmasker* help file before getting started
`./windowmasker -help`
#Review the Canada lynx genome input file (mLynCan4.pri.cur.20190110.fasta)

head mLynCan4.pri.cur.20190110.fasta
    > scaffold_11_arrow_ctg1
    > AAAAATGGCTAAGAGGGAAGGAGTAAAACATAAGAGAATCCTAAAAACTGCTGAAATGAA
    > ACACTATTCCTGAACTACACTGTAACATGAGTCCTCAACACAGCCCTCGGTGGAAATCAC
    > TACACTGCTGGAGATTCAGTCACTGCCACTGCCTTATCCCAAAGCTTACACTACACCTAA
    > CCCTAACCCTAACCCAAACCCTAACACTAACACTTACCCTAACCAAAGCCCTAGCTGTAG
    > CCCTAGCCCTGACCCTAACTCCAAGCCCCTACCCTCAAACTCTAAACCTAAGTCTAACCT
    > AAATCTACCCCCTCACTGCTAACCCTAATCCTAACTGTAACCCCTAAGACTTAACCCTAA
    > CCCTAACCAAAACCTAACCTTGTCCAGAGAAGGAGAAAGGAGAATATAAAGTCAATGTGT
    > AATAGATGCAGATTTACACTCAAAATTACAATTTCTGCAGGCGTATTGTGTTGAGATTTC
    > TACAGTAATCTGAATATCCTTAATGCTGGACAACTGTATACTTAAAATCTTTTAAGTGAT

#Apply *windowmasker* to the Canada lynx genome
`./windowmasker -mk_counts -in mLynCan4.pri.cur.20190110.fasta > genome.wmstat`
#getting an error that the first line of data is 100% ambiguous. review?
`./windowmasker -ustat genome.wmstat -outfmt fasta -in mLynCan4.pri.cur.20190110.fasta > genome-wm.fa`
#./windowmasker outputs a copy of the genome (genome-wm.fa) with interspersed repeats in lowercase actg

#Review the output:

head genome-wm.fa
    > scaffold_11_arrow_ctg1
    > AAAAATGGCTAAGAGGGAAGGAGTAAAACATAAGAGAATCCTAAAAACTGCTGAAATGAA
    > ACACTATTCCTGAACTACACTGTAACATGAGTCCTCAACACAGCCCTCGGTGGAAATCAC
    > TACACTGCTGGAGATTCAGTCACTGCCACTGCCTTATCCCAAAGCTTACACTacacctaa
    > ccctaaccctaacccaaacCCTAACACTAACACTTACCCTAACCAAAGCCCTAGCTGTAG
    > CCCTAGCCCTGACCCTAACTCCAAGCCCCTACCCTCAAACTCTAAACCTAAGTCTAACCT
    > AAATCTACCCCCTCACTGCTAACCCTAATCCTAACTGTAACCCCTAAGacttaaccctaa
    > ccctaaccaaAACCTAACCTTgtccagagaaggagaaaggagaatatAAAGTCAATGTGT
    > AATAGATGCAGATTTACACTCAAAATTACAATTTCTGCAGGCGTATTGTGTTGAGATTTC
    > TACAGTAATCTGAATATCCTTAATGCTGGacaactgtatacttaaaatcttttaagtgAT

#Now that windowmasker has done its job, we can index the genome

### Part 2: Indexing with *lastdb*
#Install LAST-979; go into the top-level LAST directory.

cd /Users/tanyalama/last-979
    #copy genome-wm.fa over to last-979 folder from                
    /Users/tanyalama/ncbi-blast-2.8.1+-src/c++/ReleaseMT/bin/genome-
    wm.fa and name as genome-wm2.fa

#Compile the programs:
`make`
#Copy the programs and scripts to a standard "bin" directory (using "sudo" to request administrator permissions):
`sudo make install`
#[Index](http://last.cbrc.jp/doc/lastdb.html) the genome using the command *lastdb*, making it suitable for comparison to highly similar sequences and alignment using *lastal*

#The general command format for *lastdb* is "lastdb mydb reference.fasta". This will create several files with names starting in "mydb".

#Consult the *lastdb* help file before getting started

lastdb --help 
lastdb -P8 -uNEAR -R11 -c clynx genome-wm2.fa 
    #-P8 Parallel threads, -uNEAR seeding scheme, -R11  
    preserves lowercase in the input, and additionally converts simple 
    sequence to lowercase, -c "masks" lowercase
    #lastdb will create several files with names starting with "clynx".
    #using genome-wm2.fa because we want to use the masked genome

### Part 3: Substitution and gap frequencies using *last-train*
#[last-train](http://last.cbrc.jp/doc/last-train.html) finds the rates (probabilities) of insertion, deletion, and substitutions between two sets of sequences. It thereby finds suitable substitution and gap scores for aligning them. It probabilistically aligns the sequences using some initial score parameters, then estimates better score parameters based on the alignments, and repeats this procedure until the parameters stop changing.

#The standard usage looks like this:
"last-train mydb queries.fasta" #should the query fasta be masked too?
#We'll try this with our Canada lynx genome and two distantly related felid genomes.

#This is the Florida panther assembly (Pco_k61_genomic.fna.gz) from NCBI
last-train -v -P0 --revsym --matsym --gapsym -E0.05 -C2 clynx Pco_k61_genomic.fna.gz > clynx-Pconcolor.mat
        #lynx-puma parameters: clynx-Pconcolor.mat
        #-v means verbose, show details on intermediate steps ran at 2:00PM
    
last-train -P0 --revsym --matsym --gapsym -E0.05 -C2 clynx PanTig1.fa > clynx-PanTig1.mat  
    #lynx-tiger parameters: clynx-PanTig1.mat
    
last-train -P0 --revsym --matsym --gapsym -E0.05 -C2 clynx PanOnc1.fa > clynx-PanOnc1.mat
    #lynx-jaguar parameters: clynx-PanOnc1.mat

### Part 4: Many-to-one felid-to-lynx alignments
#*[lastal](http://last.cbrc.jp/doc/lastal.html)* finds local alignments between query sequences, and reference sequences that have been prepared using lastdb.

#You can use it like this: "lastal humanDb dna*.fasta > myalns.maf" using .fa or gzip (.gz) compressed query files

#The *lastal* command reads files called dna*.fasta, compares them to "mydb", and writes alignments to a file called myalns.maf.

#Began run 2:30PM
lastal -v clynx Pco_k61_genomic.fa.gz > clynx-Pconcolor.maf 
    #try this for a single genome
    > lastal: batch size=8192
    > lastal: calculating matrix probabilities...
    > lastal: matrix lambda=0.229011
    > lastal: matrix letter frequencies (upper=reference, lower=query):
    > A   C   G   T  
    > 25  25  25  25 
    > 25  25  25  25 
    > lastal: getting E-value parameters...
    > 0.22852552944644766, 0.43337784178578154,
    > 0.17722035780026291, -0.1220983494599942,
    > 0.17722035780026291, -0.1220983494599942,
    > 0.044639007347921783, -0.44696673676441229,
    > 0.044639007347921783, -0.44696673676441229,
    > 0.044598546568475124, -0.44453908999761271
    > lastal: reading clynx...
    > lastal: reading Pco_k61_genomic.fa.gz...

#Review the output .maf and .mat files:

head clynx-Pconcolor.maf
    > # LAST version 979
    > #
    > # a=21 b=9 A=21 B=9 e=155 d=114 x=154 y=44 z=154 D=1e+06 E=207.799
    > # R=11 u=2 s=2 S=0 M=0 T=0 m=10 l=1 n=10 k=1 w=1000 t=4.36661 j=3 Q=0
    > # clynx
    > # Reference sequences=67 normal letters=2406165842
    > # lambda=0.228526 K=0.433378
    > #
    > #     A   C   G   T   M   S   K   W   R   Y   B   D   H   V
    > # A   6 -18 -18 -18   3 -18 -18   3   3 -18 -18   1   1   1

head clynx-Pconcolor.mat
    > # lastal version: 979
    > # maximum percent identity: 100
    > # scale of score parameters: 4.5512
    > # scale used while training: 91.024
    > 
    > # lastal -j7 -E0.05 -S1 -C2 -P0 -v -r5 -q5 -a15 -b3 clynx /var/folders/d3/4f6lh5bn3k3_2xw9p_jz7xgh0000gn/T/tmpEE4UoF
    > 
    > # aligned letter pairs: 889161.3
    > # deletes: 6499.6474
    > # inserts: 6913.2512

#Visualize your alignment using *last-dotplot*

last-dotplot clynx-Pconcolor.maf clynxPconcolor.png
    #missing PIL package, need to install

#For later: this format is for *many-to-one* genome alignment, which we will use for our ~50 lc whole genomes
#lastal -m50 -E0.05 -C2 -p hg38-panTro5.mat hg38-NEAR panTro5.fa | last-split -m1 > hg38-panTro5-1.maf
#will result in a clynx-reference.maf.gz file

#*lastal* was the slowest process in this script. We can "easily" parallelize it by processing each query sequence separately (in parallel). But each process uses quite a lot of memory, so take care that multiple parallel runs don't exceed your memory.

### Part 4: Check quality of resultant .maf files in Galaxy
#View data, E=3.5e-186, which to my understanding means that the alignment is significant (unlikely to have occurred by chance). Still unsure how to interpret EG2, lambda and K...
#lowercase regions are "masked"

    > # LAST version 979
    > #
    > # a=21 b=9 A=21 B=9 e=155 d=114 x=154 y=44 z=154 D=1e+06 E=207.799
    > # R=11 u=2 s=2 S=0 M=0 T=0 m=10 l=1 n=10 k=1 w=1000 t=4.36661 j=3 Q=0
    > # clynx
    > # Reference sequences=67 normal letters=2406165842
    > # lambda=0.228526 K=0.433378
    > #
    > #     A   C   G   T   M   S   K   W   R   Y   B   D   H   V
    > # A   6 -18 -18 -18   3 -18 -18   3   3 -18 -18   1   1   1
    > # C -18   6 -18 -18   3   3 -18 -18 -18   3   1 -18   1   1
    > # G -18 -18   6 -18 -18   3   3 -18   3 -18   1   1 -18   1
    > # T -18 -18 -18   6 -18 -18   3   3 -18   3   1   1   1 -18
    > # M   3   3 -18 -18   3   0 -18   0   0   0  -2  -2   1   1
    > # S -18   3   3 -18   0   3   0 -18   0   0   1  -2  -2   1
    > # K -18 -18   3   3 -18   0   3   0   0   0   1   1  -2  -2
    > # W   3 -18 -18   3   0 -18   0   3   0   0  -2   1   1  -2
    > # R   3 -18   3 -18   0   0   0   0   3 -18  -2   1  -2   1
    > # Y -18   3 -18   3   0   0   0   0 -18   3   1  -2   1  -2
    > # B -18   1   1   1  -2   1   1  -2  -2   1   1  -1  -1  -1
    > # D   1 -18   1   1  -2  -2   1   1   1  -2  -1   1  -1  -1
    > # H   1   1 -18   1   1  -2  -2   1  -2   1  -1  -1   1  -1
    > # V   1   1   1 -18   1   1  -2  -2   1  -2  -1  -1  -1   1
    > #
    > # Coordinates are 0-based.  For - strand matches, coordinates
    > # in the reverse complement of the 2nd sequence are used.
    > #
    > # name start alnSize strand seqSize alignment
    > #
    > # batch 0
    > a score=1986 EG2=3.4e-180 E=3.5e-186
    > s Super_Scaffold_2 137879622 363 + 146106016 ggggagcGCTTCTGCTGACTTCAGTTCACTAATTCCATTTGGTCATACTCGTCACCTTGTCTGTGTTAAGAAGGGGGCTGAGGGAGCTTGAGATGAATAGGTCGTATGCCTGCCTTCAAGAAGAGTTCCCTGGTAGTGGGAGCCTGAGCTGCAGTTAGGTGAGCCTCTGGGTGCGTCCCCTTACAGGGCAAAGAGACTCTGCAGGTGTGTTTAAGGTCAAGGGCCCTGCGATCAGAGACGATCCTGGCTTATCCAGCGGAGCCCAGTGTCATCACTGGGGTCCTCACACGAGGGAAgccggagggtcagagagagtttGGAAGGACCCCACGCTGCAGGCTTGGGAGACAGAGGACGGGGCCA
    > s PSOM01000001.1           2 363 +       564 GGGGGGCGCTTCTGCTGACTTCAGTTCACTAATTCCATTTGGTCATACTCGTCACCTTGTCTGTGTTAAGAAGGGGGCTGAGGGAGCTTGAGATGAATAGGTCGTATGCCTGCCTTCAAGAAGAGTTCCCTGGTAGTGGGAGCCCGAGCTGCAGTTAGGTGAGCCTCTGGGTGCGTCCCCTTGCAGGGCAAAGAGACTCTGCAGGGGTGTTTAAGGTCAAGGGCCCTGCAATCAGAGACGATCCTGGCTTATCCAGCGGAGCCCAGTGTGATCACTGGGGTCCTCACACGAGGGAAgccggagggtcagagagagtttGGAAGGACCCCACGCTGCAGGCTTGGGAGACAGAAGAAGGGGCCA
    > 
    > a score=690 EG2=1.4e-51 E=3e-57
    > s Super_Scaffold_2 137879985 131 + 146106016 TGTCCCAGACGTCCGACCTCCAGGCCTCCAAGCTACTGAGCGTGCGGCAACTCGTTAAGGCCGCCGTGGGAAGCTGACACAGGCACAGTAAATGTGTGCAGGGGGCCAGCGCGCAGGAAGGGGTGCAGggt
    > s PSOM01000001.1         433 131 +       564 TGTCCCAGACGTCCGACCTCCAGGCCTCTGAGCTACTGAGCGTGCGGCAACTCGTTAAGGCCGCCGTGGGAAGCTGACACAGGCACAGGAAATGCGTGCAGGGGGCCAGCGCGCAGGAAGGGGTGCAGGGT
    > 
    > a score=156 EG2=1.4e+02 E=0.00037
    > s Super_Scaffold_12 157953951 26 + 158927374 GGCCCGTGGCCCCTTCTTCTGTCTCC
    > s PSOM01000001.1          193 26 -       564 GGCCCGTGGCCCCTTCTTCTGTCTCC
    > 
    > a score=1419 EG2=6.4e-124 E=8.2e-130
    > s Super_Scaffold_4 42856003 282 + 220954660 TAAAGCAGTATTTGTCATGTGTTTACAGTTACAATATGGACTGAATTTATAGTGATGCCAtgacataaaatgaaaagatctCTTGAAAGATAAATTCCTGATATCTGCCTACCATTCTCATAGAAGCAGAATTCAGGGGGAAAGTCTGCGTTTTG----GTAAGgatttactgttttgttttgttttgtttggataaAtaagaatttactgttttAAGAAACTAGACAATGAATATACatgtaaaggaaaattaaatttaattttctttgatcctttcttcttttt
    > s PSOM01000002.1         55 281 +       520 taaagcagtatttGTCATGTGTTTACAGTTACAATATGGACTGAATTTATACTGATGCCAtgacataaaatgaaaagatctcTTGAAAGATAAATTCCTGATATCTGCCTGCCATTCTCATAGAAGCAGAATTCAGGGGGAAAGTCTGTGTTTTGGTAAGTAAGgatttac-----tgttttgttttgtttggataaataagaatttactgttttAAGAAACTAGACAATGAATATGCATGTAAAGgaaacttaaatttaattttctttgatcctttcttcttttt
    > 
    > a score=1152 EG2=2e-97 E=3.1e-103
    > s Super_Scaffold_4 42856512 200 + 220954660 atcctttcttctttttgactCTCCTTTTATTCAGAAACAGTGCCTTGCATATATTGGCATTAAATTCATGCAtgttgtgggtgcctgggtggctcagttggttaagcaagcatctgacttcggctcaggtcatgatcgcttggtttgtgagtttgagccctgcgttgggctctgtgctgacagctcagagcctggagcct
    > s PSOM01000002.1        320 200 +       520 atcctttcttctttttgactCTCCTTTTATTCGGAAACAGTGCCTTGCATATATTGGCATTAAATTCATGCAtgttgtgggtgcctgggtggctcagttggttaagcaagcatttgacttcggctcaggtcatgatcgcttggtttgtgagtttgagccctgcgttgggctctgtgctgacagctcagagcctggagcct
    > 
    > a score=6411 EG2=0 E=0
    > s Super_Scaffold_14 109755312 1185 + 115715915 tctcaaaaatgaataaacatttaagaaaaaaaaattaaaaaacccaaaaactcaaGAGTTGAAAAGAATTTTCTGAACTGATAAAGAatgtttcaaaaaatgaaaatagctgtCAACTGTTGAAACACTACAACActccttttaaaatcagaaagaagaccaaaataaaataaaataaaaattagaaagaagacaaaaatgctCCCTATTTCCACTTCTGATCAACAATTTACTGAAGGTCTTAGGTTTTGGAAGGAAAcgagaagaaaaatgtaaagtttaggcattcaaaaggaagaaacaaaactctt--TGTTGCAAATGATATGACTAGGCAAGCCAAACAAATCTACAAATTGTTGAAATTAATAAAGAGTTTCACAGGGTTGAATAAAGTCCCAAGAGGGTTTTTCATGGAACTTGACAAGCTGATTGCAAAATTTAGATGGAAGCATAAATGATTATGATGAAGAGCCAAGGCAGTCCTGAAGCAAAACCGGTTGGGGGACTCATTCCATCAGATACCAAAATTTATTCTAAGATTATTAAGTTCATTATAGCCCTAGGGCAGGAATAGACAAACTGACCAGTGGAATGAACCAAGAACCAGGCATTGACTCTCGTAGATGATAGAAGGGACACGGTAGGCCACTGGAGAAGGGAGGGACCATTTAATAAGTCATGCTGAGACGACAGGTTCTCCTTATGGAGAGCCAAACTGGATTTCTACCTTCTGCAATACATAGTCAACACCACGAGATTAAAGTCTTACATAGGAAAGGCAAAATTTTAAACTGGGAgaagaatatctttatgacctcATAGTGAAGGAGAATTTCTCTAAGCACAAGAAGCAGCAGCTATCAGGGAGAAGATCGGTCCTTTTAACTACATTCAAAATagctttgggatgcctgggtggcccagttggttaagcgtctgactcttggtttcggctcaagtcgtgatctcacagttagtgggtcaggctccacacagatagtgcagagcctgcttgggattctctctctccctctctctctgcccctccccacctgaaaataaataaataaaaactttaaaaagtaattggcTTTTATTGGAGaacctggctggctaagtcagtgGAGGGTGAGGCTGTTGATCTCAGTATTGTAaattcaagcaccatgttgggtgtagagatt
    > s PSOM01000003.1           36 1181 +      1227 tctcaaaaatgaataaacatttaagaaaaaaaaattaaaaaacccaaaaactcaagagttgaaaagaattttcttaactgataaagaatatttcaaaaaatgaaaatagctgtCAACCGTTGAAACACTACAACActccttttaaaatcagaaagaagaccaaaataaagtaaaataaaaattagaaagaagacaaaaatgctCCCTATTTCCACTTCTGATCAACAATTTACTGAAGGTCTTAGGTTTTGGAAGGAAAcgagaagaaaaatgtaaagtttaggcattcaaaaggaagaaacaaaactcttCGTGTTGCAAATGATATGACTAGGCAAGGCAAACAAATCTACAAATTGTTGAAATTAATAAAGAGTTTCACAGGGTTGAATAAAGTCCCAAGAGGGTTTTTCATGGAACTTGACAAGCTGATTGCAAAACTTAGATGGAAGCATAAATGATCATGATGAAGAGCCAAGGCAGTCCTGAAGCAAAACCGGTTGGGGGATTCACTCCATCAGATACCAAAATTTATTCTAAGATTATTAAGGTCAGTATAGTCCTGGGGCAGGAATAGACAAACTGACCAGTGGAATGAACCAAGAACCAGGCATTGACTCTCGTAGATGATAGAAGGGACATGGTAGGCCACTGGAGAAGGGAGGGACCATTTAATAAGTCATGCTGAGACAACAGGTTCTCCTTATGGAGAGCCAAACTGGATTTCTACCTTCCGCAATACATAGTCAACACCACGAGATTAAAGTCTTACATAGGAAAGGCAAAATTTTAAACTGGGAGAAGAATA----TATGACCTCATAGTGAAGGAGAATTTCTCTAAGCACAAGAAGCACCAGCTATCAGGGAAAAGATCAGTCCCTTTAACTACATTCAAAATagctttgggatgcctgggtggctcagttggttaagcgtctgactcttggtttcggctcaagtcgtgatctcacagttagtgggtcaggctccacacagatagtgcagagcctgcttgggattctctctctccctctctctctgcccctccccacctgaaaataaataaataaaaactttaaaaagtaattggcTTTCATTGGAGAACCTGGCTGGCTAAGT--GTGGAGGGTGAGGCTGTTGATCTCAGTATTGTAaattcaagcaccatgttgggtgtagagatt
    > 
    > 

#Note: it's occurred to me that we masked the reference sequence (lynx) and mapped the other felids *to* lynx. We'll want to re-do this by mapping the lynx *to* one of the well-annotated reference genomes (probably *Felis catus*). Basically, lynx should be the query fasta and *Felis catus* should be the reference fasta

#Next steps: we actually need a bam file for some downstream analyses (including PSMC) so we need to re-do this and convert .maf to sam/bam...?