# Scripts for the Publication:

Ullrich KK, Tautz D

## Data sources

Genome mapping files for _Mus musculus domesticus GER_, _Mus musculus domesticus FRA_, _Mus musculus domesticus IRA_, _Mus musculus musculus AFG_, _Mus musculus castaneus CAS_ and _Mus spretus SPRE_ were obtained from <http://wwwuser.gwdg.de/~evolbio/evolgen/wildmouse/m_m_domesticus/genomes_bam/>, <http://wwwuser.gwdg.de/~evolbio/evolgen/wildmouse/m_m_musculus/genomes_bam/>, <http://wwwuser.gwdg.de/~evolbio/evolgen/wildmouse/m_m_castaneus/genomes_bam/>, <http://wwwuser.gwdg.de/~evolbio/evolgen/wildmouse/m_spretus/genomes_bam/>.

For mapping details please look into the original publication ([Harr et al. 2016](http://www.nature.com/articles/sdata201675)) <http://www.nature.com/article-assets/npg/sdata/2016/sdata201675/extref/sdata201675-s7.docx>.

## Get masking regions for individual samples and natural populations

For masking genomic regions in natural populations which showed low coverage based on the genomic mapping BAM files we only considered the stable chromosomes from the reference GRCm38 _mm10_ <http://www.ncbi.nlm.nih.gov/projects/genome/assembly/grc/mouse/>.

The BAM files were processed with 'genomeCoverageBed' to obtain site specific genome coverage and further united with 'unionBedGraphs'. The per population combined coverage was further processed to only retain regions with a coverage smaller than 5 resulting as the masking regions.

genomeCoverageBed example for the population _Mus musculus musculus AFG_:
```
#example for populaton Mmm_AFG:
#
#used BAM files:
#
#AFG1_396.bam
#AFG2_413.bam
#AFG3_416.bam
#AFG4_424.bam
#AFG5_435.bam
#AFG6_444.bam

$REFERENCE=mm10.fasta

for file in *.bam; do genomeCoverageBed -ibam $file -bga -g $REFFERENCE > $file".bga";done
```

unionBedGraphs example for the population _Mus musculus musculus AFG_:
```
INPUT1=AFG1_396.bam.bga
INPUT2=AFG2_413.bam.bga
INPUT3=AFG3_416.bam.bga
INPUT4=AFG4_424.bam.bga
INPUT5=AFG5_435.bam.bga
INPUT6=AFG6_444.bam.bga

OUTPUT=Mmm_AFG.combined.bga

unionBedGraphs -i $INPUT1 $INPUT2 $INPUT3 $INPUT4 $INPUT5 $INPUT6 | awk -v OFS='\t' 'BEGIN {sum=0} {for (i=4: i<=NF; i++) sum+=$1; print $1,$2,$3,sum; sum=0}' > $OUTPUT
```

get masking region example for the population _Mus musculus musculus AFG_:
```
INPUT=Mmm_AFG.combined.bga
OUTPUT=Mmm_AFG.combined.bga.stcov5
awk '{if($4<5) print $0}' $INPUT > $INPUT".stcov5"
bedtools merge -i $INPUT".stcov5" > $INPUT".stcov5.merge"
awk -v OFS='\t' '{print $1,$2,$3,4}' $INPUT".stcov5.merge" > $OUTPUT
```

_used software:_

+ bedtools v2.24.0 <http://bedtools.readthedocs.io/en/latest/>
+ awk

## SNP and INDEL calling

For SNP and INDEL calling the BAM files were processed with 'samtools mpileup' and 'bcftools call' with relaxed quality options to retain information in CNV regions.
Due to the amount of data first mpileup files were generated for each BAM file and merged with 'bcftools merge'. After merging mpileup files, 'bcftools call' was used to generate the final VCF file.

samtools mpileup | bcftools call example for the population _Mus musculus musculus AFG_:
```
#example for population Mmm_AFG:
#
#used BAM files:
#
#AFG1_396.bam
#AFG2_413.bam
#AFG3_416.bam
#AFG4_424.bam
#AFG5_435.bam
#AFG6_444.bam

#generating indiviual mpileup files for each BAM file

$REFERENCE=mm10.fasta

INPUT1=AFG1_396.bam
INPUT2=AFG2_413.bam
INPUT3=AFG3_416.bam
INPUT4=AFG4_424.bam
INPUT5=AFG5_435.bam
INPUT6=AFG6_444.bam

samtools mpileup -q 0 -Q 10 -A -d 99999 -t DP,AD,ADF,ADR -vf $REFERENCE -o $INPUT1".mpileup.q0Q10.vcf" $INPUT1
samtools mpileup -q 0 -Q 10 -A -d 99999 -t DP,AD,ADF,ADR -vf $REFERENCE -o $INPUT2".mpileup.q0Q10.vcf" $INPUT2
samtools mpileup -q 0 -Q 10 -A -d 99999 -t DP,AD,ADF,ADR -vf $REFERENCE -o $INPUT3".mpileup.q0Q10.vcf" $INPUT3
samtools mpileup -q 0 -Q 10 -A -d 99999 -t DP,AD,ADF,ADR -vf $REFERENCE -o $INPUT4".mpileup.q0Q10.vcf" $INPUT4
samtools mpileup -q 0 -Q 10 -A -d 99999 -t DP,AD,ADF,ADR -vf $REFERENCE -o $INPUT5".mpileup.q0Q10.vcf" $INPUT5
samtools mpileup -q 0 -Q 10 -A -d 99999 -t DP,AD,ADF,ADR -vf $REFERENCE -o $INPUT6".mpileup.q0Q10.vcf" $INPUT6

bgzip $INPUT1".mpileup.q0Q10.vcf"
bgzip $INPUT2".mpileup.q0Q10.vcf"
bgzip $INPUT3".mpileup.q0Q10.vcf"
bgzip $INPUT4".mpileup.q0Q10.vcf"
bgzip $INPUT5".mpileup.q0Q10.vcf"
bgzip $INPUT6".mpileup.q0Q10.vcf"

tabix $INPUT1".mpileup.q0Q10.vcf.gz"
tabix $INPUT2".mpileup.q0Q10.vcf.gz"
tabix $INPUT3".mpileup.q0Q10.vcf.gz"
tabix $INPUT4".mpileup.q0Q10.vcf.gz"
tabix $INPUT5".mpileup.q0Q10.vcf.gz"
tabix $INPUT6".mpileup.q0Q10.vcf.gz"

#merge individual mpileup files

MPILEUPLIST=AFG.mpileup.list

echo $INPUT1".mpileup.q0Q10.vcf.gz" >> $MPILEUPLIST
echo $INPUT2".mpileup.q0Q10.vcf.gz" >> $MPILEUPLIST
echo $INPUT3".mpileup.q0Q10.vcf.gz" >> $MPILEUPLIST
echo $INPUT4".mpileup.q0Q10.vcf.gz" >> $MPILEUPLIST
echo $INPUT5".mpileup.q0Q10.vcf.gz" >> $MPILEUPLIST
echo $INPUT6".mpileup.q0Q10.vcf.gz" >> $MPILEUPLIST

MPILEUPOUTPUT=AFG.mpileup.q0Q10.vcf.gz

bcftools merge -m all -O z -o $MPILEUPOUTPUT -l $MPILEUPLIST

#call SNP and INDEL

OUTPUT=AFG.mpileup.q0Q10.bcfcall.mv.vcf.gz

bcftools call -O z -f GQ -m -v -o $OUTPUT $MPILEUPOUTPUT
```

_used software:_

+ samtools 1.3.1-36-g613501f (using htslib 1.3.1-59-g0f2a88a)
+ bcftools 1.3.1-39-gd797e86 (using htslib 1.3.1-59-g0f2a88a)
+ bgzip v1.3
+ tabix v1.3

## K80 distance calculation

### Get population specific SNPs

_used software:_

### Calculate population specific Consensus sequence

_used software:_

### Calculate K80 distance between populations

## Dxy distance calculation

### Calculate Dxy distance between populations

_used software:_

### Calculate Dxy distance between individuals and populations

_used software:_