---
# System prepended metadata

title: Mapping ddRAD reads to a reference genome
tags: [RADSeq]

---

---
title: Mapping ddRAD reads to a reference genome
breaks: false
tags: RADseq, Nora
---


# Mapping ddRAD reads to a reference genome
### Preparing the reference:
Good idea to make a new directory for your reference files. The reference needs to be uncompressed for this analysis. Unzip the reference:
```
mkdir ./reference
mv CAJRAG01.1.fsa_nt.gz ./reference/
gunzip CAJRAG01.1.fsa_nt.gz
```
Index the reference genome:
```
bwa index CAJRAG01.1.fsa_nt
```
This will create a few other files with the same basename. They all need to be in the same directory with the .fsa_nt file in order for mapping and such to work.

### Map the reads:
Run from the directory with the read files, this loop will align the reads for each pair of read files, generating a .sam file
```
bwa mem -M -t 4 /home/bioinfo8/reference_mapping/reference_files/CAJRAG01.1.fsa_nt sampleID.1.fq.gz sampleID.2.fq.gz > sampleID.sam
```

For references that might be more distantly related, you might get low percentage of reads mapping to the reference. You can check the percent mapping to the reference with:
```
samtools flagstat sampleID.sam
```
Good mapping percentages are in the 90's. You can mess around with the -K values of bwa mem to see if you get a higher percent mapping to the reference. This is a tradeoff because the lower the -K, that means that the fewer exact basepairs need to match the reference in order to start mapping a read. Default is 19 when bwa mem is run without a -K value. It can be a good idea to run bwa mem on one sample with several values of -K : default, -K 16 ; -K 14; -K 12 ; -K 10. Evaluate which maps the most reads and then use the loop for the rest of the samples:
```
for i in *.1.fq.gz; do
NM=${i//.1.fq.gz/}
REVERSE=${i//.1.fq.gz/.2.fq.gz}
bwa mem -M -t 4 -K 10 /home/bioinfo8/reference_mapping/reference_files/CAJRAG01.1.fsa_nt $i $REVERSE > $NM.sam ; done
```
Now convert the .sam files into .bam files:
```
for i in *.sam; do
NM=${i//.sam/}; 
samtools view $i -b -o $NM.bam; done
```
and then sort the .bam files:
```
for i in *.bam; do  NM=${i//.bam/}; samtools sort $i -o $NM.sort.bam; done
```
As long as your sorted .bam files are looking good, then you can delete the .sam files because they are quite large and take up a lot of space. I put all my sorted .bams into their own folder.

### Build stacks loci from aligned reads:
Good idea to do this in screen
```
gstacks -I ~/sorted_bams/ -M ~/popmap -O ~/reference_stacks/ -t 8
```
#note from May 2024: maybe should not use --rm-pcr-ducplicates https://groups.google.com/g/stacks-users/c/qmAxcH7cbus

and now run populations, also in screen
```
populations -P ~/reference_stacks/ -M ~/popmap -r 0.65 --vcf --genepop --fstats --smooth --ordered-export --hwe -t 8
```
This is where I stopped as far as looking for amplicon loci with the reference mapped stacks, because there were many SNPS from the same loci. From what I can gather, the flag --write-random-snp doesn't work with reference aligned stacks, and it is suggested instead to used --ordered-export (also note that this is --ordered_export in the STACKS documentation, but it runs only with --ordered-export with "-" instead of "_"  go figure)

## Summarizing results:
To calculate average missingness in the filtered vcf
```
vcftools --gzvcf engleri_filtered1.recode-1.vcf --missing-indv --out dp
vcftools --gzvcf engleri_filtered1.recode-1.vcf --depth --out dp
```
