test_saf_sfs

Below is the popARRAY script for Togiak and Kotzebue pairwise Fst in a Manhattan plot

# togiak saf
angsd -b /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/pherr_togiak_bamslist.txt \
-r ${chrom}: \
-sites /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/gls/togiak/pherr_wholegenome_polymorphic.sites \
-ref /center1/GLASSLAB/salmgren/atlantic_herring/GCF_900700415.2_Ch_v2.0.2_genomic.fna \
-anc /center1/GLASSLAB/salmgren/atlantic_herring/GCF_900700415.2_Ch_v2.0.2_genomic.fna \
-out /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/diversity/pherr_${chrom}_TG_polymorphic_folded \
-nThreads 10 \
-GL 1 \
-doGlf 3 \
-doMaf 1 \
-doMajorMinor 1 \
-trim 0 \
-C 50 \
-minMapQ 15 \
-minQ 15 \
-doCounts 1 \
-setminDepth 20 \
-setmaxDepth 100 \
-remove_bads 1 \
-uniqueOnly 1 \
-only_proper_pairs 1 \
-doSaf 1 

# kotzebue saf
angsd -b /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/pherr_kotzebue_bamslist.txt \
-r ${chrom}: \
-sites /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/gls/kotzebue/pherr_wholegenome_polymorphic.sites \
-ref /center1/GLASSLAB/salmgren/atlantic_herring/GCF_900700415.2_Ch_v2.0.2_genomic.fna \
-anc /center1/GLASSLAB/salmgren/atlantic_herring/GCF_900700415.2_Ch_v2.0.2_genomic.fna \
-out /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/diversity/pherr_${chrom}_KZ_polymorphic_folded \
-nThreads 10 \
-GL 1 \
-doGlf 3 \
-doMaf 1 \
-doMajorMinor 1 \
-trim 0 \
-C 50 \
-minMapQ 15 \
-minQ 15 \
-doCounts 1 \
-setminDepth 19 \
-setmaxDepth 95 \
-remove_bads 1 \
-uniqueOnly 1 \
-only_proper_pairs 1 \
-doSaf 1 
    # -doSaf estimates sample allele frequency (saf)

# realSFS (program) uses saf files produced by angsd -doSaf
        #estiamte saf, then do sfs
realSFS /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/diversity/pherr_${chrom}_TG_polymorphic_folded.saf.idx \
-fold 1 /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/diversity/pherr_${chrom}_KZ_polymorphic_folded.saf.idx \
-fold 1 > /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/fst/pherr_${chrom}_TG-KZ_polymorphic_folded.sfs
        # 2D saf into sfs

# calculate per-site Fst
realSFS fst index /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/diversity/pherr_${chrom}_TG_polymorphic_folded.saf.idx \
-fold 1 /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/diversity/pherr_${chrom}_KZ_polymorphic_folded.saf.idx \
-fold 1 -sfs /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/fst/pherr_${chrom}_TG-KZ_polymorphic_folded.sfs \
-fstout /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/fst/pherr_${chrom}_TG-KZ_polymorphic_folded.sfs.pbs \
-whichFst 1
        # which Fst estimator 

# sliding window analysis
realSFS fst stats2 /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/fst/pherr_${chrom}_TG-KZ_polymorphic_folded.sfs.pbs.fst.idx \
-win 1 -step 1 > /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/fst/pherr_${chrom}_TG-KZ_polymorphic_folded.sfs.pbs.fst.txt
        # across every site, sliding window of 1

# concateante fst files for each chrom 

Plot from concatenated fst output files from popARRAY script (R plotting script is the same for both)

setwd("/Users/sydneyalmgren/Documents/SA_P_Herring/bering_sea_pop_gen_output/sfs_fst_spawning")

# read in chromosome list 
chroms <- as.data.frame(read.csv("chromosome_list.csv"))
chroms <- as.data.frame(chroms)
# make sure that numbered column in chromosome list (header for 1,2,3,...) is CHR for easier use with qqman

### calculate by chrom, concatenate
### not necassarily correct? saf and sfs steps not in order 
tg.kz1 <- read_delim("pherr_TG-KZ_wholegenome_fst.txt", 
                    skip = 1, 
                    col_types = cols(), 
                    col_names = c("region", "chr", "midpos", "nsites", "fst"), 
                    delim = "\t")
# change to dataframe
tg.kz1 <- as.data.frame(tg.kz1)

# filter out sites <0, make MB pos column 
tg.kz1 <- filtered_tg.kz1 <- tg.kz1 %>%
  filter(fst >= 0)
tg.kz1 <- tg.kz1 %>%
  mutate(midpos_Mb = midpos/1000000)

# add chroms to fst dataframe
tg.kz1 <- left_join(tg.kz1, chroms, by = "chr")

### plot with qqman package 
manhattan(tg.kz1, chr = "CHR", bp = "midpos", p = "fst", snp = "nsites", 
          main = "Togiak vs Kotzebue (calculated by chrom, concatenated)",
          logp = FALSE, cex = 0.5, cex.axis = 0.8, ylab=expression(italic(F)[ST]))

This is two scripts: one that estimates saf for Togaik and Kotzebue separately and calculates folded sfs, then calculates Fst. These scripts are not submitted as an array

# togiak saf 
angsd -b /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/pherr_togiak_bamslist.txt \
-sites /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/gls/togiak/pherr_wholegenome_polymorphic.sites \
-ref /center1/GLASSLAB/salmgren/atlantic_herring/GCF_900700415.2_Ch_v2.0.2_genomic.fna \
-anc /center1/GLASSLAB/salmgren/atlantic_herring/GCF_900700415.2_Ch_v2.0.2_genomic.fna \
-out /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/diversity/togiak/pherr_TG_wholegenome_polymorphic_folded \
-nThreads 10 \
-GL 1 \
-doGlf 3 \
-doMaf 1 \
-doMajorMinor 1 \
-trim 0 \
-C 50 \
-minMapQ 15 \
-minQ 15 \
-doCounts 1 \
-setminDepth 20 \
-setmaxDepth 100 \
-remove_bads 1 \
-uniqueOnly 1 \
-only_proper_pairs 1 \
-doSaf 1 

# convert to SFS
realSFS -fold 1 /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/diversity/togiak/pherr_TG_wholegenome_polymorphic_folded.saf.idx > /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/diversity/togiak/pherr_TG_wholegenome_polymorphic_folded.sfs

## repeated above for Kotzebue


# 2D SFS [pop1] [pop2] > [pop1_2]
realSFS /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/diversity/togiak/pherr_TG_wholegenome_polymorphic_folded.saf.idx \
-fold 1 /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/diversity/kotzebue/pherr_KZ_polymorphic_folded.saf.idx \
-fold 1 > /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/diversity/pherr_wholegenome_polymorphic_TG-KZ_folded.sfs

# Fst binary files
realSFS fst index /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/diversity/togiak/pherr_TG_wholegenome_polymorphic_folded.saf.idx \
-fold 1 /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/diversity/kotzebue/pherr_KZ_polymorphic_folded.saf.idx \
-fold 1 -sfs /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/diversity/pherr_wholegenome_polymorphic_TG-KZ_folded.sfs \
-fstout /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/fst/pherr_wholegenome_polymorphic_TG-KZ_folded.sfs.pbs \
-whichFst 1

# get Fst values by sliding window 
realSFS fst stats2 /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/fst/pherr_wholegenome_polymorphic_TG-KZ_folded.sfs.pbs.fst.idx \
-win 1 -step 1 > /center1/GLASSLAB/salmgren/lcwgs/ebs_popgen/fst/pherr_wholegenome_polymorphic_TG-KZ_folded.sfs.pbs.fst.txt

Plot for non-array saf, sfs, fst script

### calculate whole genome saf and sfs, then fst folded with whole genome 
tg.kz2 <- read_delim("pherr_wholegenome_polymorphic_TG-KZ_folded.sfs.pbs.fst.txt", 
                    skip = 1, 
                    col_types = cols(), 
                    col_names = c("region", "chr", "midpos", "nsites", "fst"), 
                    delim = "\t")
# change to dataframe
tg.kz2 <- as.data.frame(tg.kz2)

# filter out sites <0, make MB pos column 
tg.kz2 <- filtered_tg.kz2 <- tg.kz2 %>%
  filter(fst >= 0)
tg.kz2 <- tg.kz2 %>%
  mutate(midpos_Mb = midpos/1000000)

# add chroms to fst dataframe
tg.kz2 <- left_join(tg.kz2, chroms, by = "chr")

### plot with qqman package 
manhattan(tg.kz2, chr = "CHR", bp = "midpos", p = "fst", snp = "nsites", 
          main = "Togiak vs Kotzebue (calculated by whole genome)",
          logp = FALSE, cex = 0.5, cex.axis = 0.8, ylab=expression(italic(F)[ST]))

These Manhattan plots look very different?? Also I can’t figure out where along the way the global Fst is estiamted (I’ve done this before, but I didn’t take good enough notes)