---
# System prepended metadata

title: Andrew Kittentails Project
tags: [RADSeq, ' Andrew']

---

---
title: Andrew Kittentails Project
break: false
tags: RADSeq, Andrew
---

# Synthyris bullii sequence data
This document contains notes for a population genetics assessment for Synthyris bullii (Kittentails) using the STACKS pipeline with and without a reference genome. From across Illinois and Indiana, 120 individuals were sequenced using NovaSEQ NextGen Sequencing (10 individuals from 12 populations). The reference genome was generated from an individual at the Nachusa Doug population at Nachusa Grasslands, Illinois, was sequenced by PacBio at Mt. Sinai using a SMRTcell on Revio, and assembled using HifiASM.

–adapter enzyme 1 (P1-EcoRI): ACACTCTTTCCCTACACGACGCTCTTCCGATCT
–adapter enzyme 2 (P2-MspI): GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT
        



# ADMIXTURE

* requires low missing data
* pruned for Linkage disequilibrium

## Pipeline Summary
- Filter missing data (--max-missing 0.85)
- Convert to PLINK
- Prune for LD (--indep-pairwise)
- Run ADMIXTURE across K values
- Visualize ancestry proportions


#### 1.  begin with vcf file used for pop gen
vcf < - vcfR.gq20.smiss96.70.vcf

```
**** Object of Class vcfR *****
85 samples
163 CHROMs
2,327 variants
Object size: 11.5 Mb
18.5 percent missing data
*****        *****         *****
```

#### 2. Filter Missingness in R
```
# Keep SNPs with <20% missing
keep_snps <- which(missing_snp < 0.2)
vcf_filtered <- vcf[keep_snps, ]

# Optionally filter samples too
keep_samples <- which(missing_sample < 0.2)
vcf_filtered@gt <- vcf_filtered@gt[, c(1, keep_samples + 1)]  # +1 for FORMAT column

# save filtered vcf

write.vcf(vcf_filtered, file = "filtered.vcf")

```
#### 3. Convert VCF to PLINK
```
plink --vcf filtered.vcf --make-bed --out filtered --allow-extra-chr
```

#### 4. Replace 'chromosome' labels with integers
* Extract unique chromosomes from filtered vcf file
```
bcftools query -f '%CHROM\n' filtered.vcf | sort -u > scaffolds.txt
```
* Create integer-chromosome map 
```
awk '{print $1 "\t" NR}' scaffolds.txt > integer_to_chr_map.txt
```
* replace 'chromosomes' with integers in .bim file
```
awk 'NR==FNR{a[$1]=$2; next} {$1=a[$1]; print}' integer_to_chr_map.txt filtered.bim > filtered.updated.bim
mv filtered.bim contig.filtered.bim
mv filtered.updated.bim filtered.bim
```

#### 3. LD Pruning with PLINK moving window
* find snps in LD and write to pruned files:
    * window size = 50
    * step size (in SNPs) = 10
    * r-squared threshold = 0.2
* prune  from original plink file
```
plink --bfile filtered --indep-pairwise 50 10 0.2 --out pruned --allow-extra-chr
plink --bfile filtered --extract pruned.prune.in --make-bed --out data_pruned --allow-extra-chr
```






### If you lost a lot of data in the above, try imputation before ADMIXTURE...

#### 2A. remove rare SNPs in R (that only occur once, OPTIONAL, SKIP FIRST TIME THROUGH))
 * Jump to step X. filter missing data first time through*

```
chrom_counts <- table($vcf$@fix[, "CHROM"])
keep_chroms <- names(chrom_counts[chrom_counts > 1])
vcf_rare_snp_removed <- vcf[vcf@fix[, "CHROM"] %in% keep_chroms, ]
write.vcf(vcf_rare_snp_removed, file = "$SAVEPATH$")
```

#### 3A. Impute missing genotypes with Beagle in terminal (OPTIONAL, SKIP FIRST TIME THROUGH)
I ended up losing 33 samples and about 1500 snps after filtering for
<20% missing by samples and <20% missing by snp, so am trying impute 
```
java -jar beagle.jar gt=filtered.vcf out=imputed
```
#### 4A. Convert VCF to PLINK
```
plink --vcf imputed.vcf --make-bed --out imputed --allow-extra-chr
```
#### 5A. Rename scaffolds to integers and create map (necessary for admixture - it doesn't like ptgs or whatever your snps are called)
```
cut -f1 imputed.bim | sort | uniq | nl > scaffold_map.txt
awk 'NR==FNR{map[$2]=$1; next} {$1=map[$1]; print}' scaffold_map.txt imputed.bim > imputed.integer.bim
```
rename original imputed and integer imputed for downstream
```
mv imputed.bim imputed.contig.bim
mv imputed.integer.bim imputed.bim
```
#### 6A. Filter SNPs by Missingness
```
plink --bfile imputed --geno 0.05 --make-bed --out imputed.filtered 
```

#### 7A. LD Pruning
```
plink --bfile imputed.filtered --indep-pairwise 50 5 0.2 --out pruned
plink --bfile imputed.filtered --extract pruned.prune.in --make-bed --out data_pruned 
```
#### 8A. Run ADMIXTURE across K values
```
for K in {4..12}; do admixture --cv pruned.bed $K | tee log${K}.out; done
```
#### 9 Make pop map based on samples from .fam file
* extract column from .fam file on command line
```
cut -f2 data_pruned.fam > sample_names.txt
```
* open .txt file, copy and paste into excel
* use =left($a$1, "number of characters in sample name")
* then = left($b$1$, "number of characters in pop code")
* repeat for all samples
* copy and paste column 2 and 3 as values only and delete other columns
* save as .csv
* or at least this worked for me with samples like HH1234. I just wanted the sample name and a column with the first 2 letters (HH)

#### ADMIXTURE PLOTS in R
```


```



# Revisiting sequence data in 2024
New things to try: 
* more thoroughly assess raw read quality with fastqc 
* use 'cut adapt' to remove excess C's generated duting 'dark cycle' of NovaSeq run prior to process rad tags
* explore 'runs of homozygosity' as an alternative to Fis as a measure of inbreeding
* do not trim reads to 100 instead of 150
* map raw reads to reference genome

### FastQC results of raw data

* For all sublibraries, the first  several reads have lower per base sequence quality and extreme (100/0%) per base sequence content between reads 6 and 10. 
    * Trimming the first 10 bps to remove these reads with cutadapt (installed with pip on quest, I don't have an allocation so can't submit jobs)
        *      cutadapt -u 10 -o /projects/p31922/01_raw_reads/trimmed_P1_S10_R1_001.fastq.gz -i /projects/p31922/01_raw_reads/P1_S10_R1_001.fastq.gz
        * running cutadapt on one core on one set of reads (73mil reads) took 10 minutes

* Trimmed reads look slightly better (at least improved Quality for first reads and better per base sequence) Running trimmomatic with a sliding scale to remove reads that have a score <20 in a 4 base window, maintaining reads over 100bp. Code:

        *    java -jar trimmomatic-0.39.jar SE -phred33 /projects/p31922/01_raw_reads/trimmed_P1_S10_R1_001.fastq.gz /projects/p31922/01_raw_reads/qualitytrimmed_P1_S10_R1_001.fastq.gz SLIDINGWINDOW:4:20 MINLEN:100 

* Not much improvement outside of quality scores - going to move forward with both datasets and compare SNP output

Process Radtags

Multiplexed Reads Notes (P1 only - can't find log file for P2/P3):
* Original raw reads truncated to 100: 137492637
* Re-multiplexed raw reads with no truncation: 137567342 (93.9%)
* Trimmed reads: 3599550 (only 2%; must've trimmed barcode - not a great idea)
* Trimmed and quality trimmed reads:  2289410 (2.1% - not a good way to go either)
* Truncating during process_radtags is perhaps the better way to deal with issues at the beginning/end of sequences

# Reference aligning sequence data onto assembled genome

Genome details: 
        
        (BUSCO v5.5; embryophyta lineage (50 genomes and 1614 BUSCOs), gene predictor: metaeuk):
        ***** Results: *****
        C:98.5%[S:74.4%,D:24.1%],F:0.4%,M:1.1%,n:1614
        1590    Complete BUSCOs (C)
        1201    Complete and single-copy BUSCOs (S)
        389     Complete and duplicated BUSCOs (D)
        7       Fragmented BUSCOs (F)
        17      Missing BUSCOs (M)
        1614    Total BUSCO groups searched

        Assembly Statistics:
        1065    Number of scaffolds
        1065    Number of contigs
        1337213455      Total length
        0.000%  Percent gaps
        9 MB    Scaffold N50
        9 MB    Contigs N50
        
        (BUSCO v5.5; eudicot lineage (31 genomes and 2326 BUSCOs), gene predictor: metaeuk):
        ***** Results: *****
        C:96.5%[S:71.4%,D:25.1%],F:0.4%,M:3.1%,n:2326
        2244    Complete BUSCOs (C)
        1661    Complete and single-copy BUSCOs (S)
        583     Complete and duplicated BUSCOs (D)
        10      Fragmented BUSCOs (F)
        72      Missing BUSCOs (M)
        2326    Total BUSCO groups searched

        Assembly Statistics:
        1065    Number of scaffolds
        1065    Number of contigs
        1337213455      Total length
        0.000%  Percent gaps
        9 MB    Scaffold N50
        9 MB    Contigs N50
   
 Index assembled genome (command line in quest took about 20 minutes):
 
     module load bwa
     bwa index asm.np2.p_ctg.fa
     
     
Mapping reads to reference with bwa mem (should i try minimap2?)

    # Define input and output directories
    input_dir="/projects/p31922/02_multiplexed_reads/1.multiplexed_raw_trunc100/realign_reads/"
    output_dir="/projects/p31922/06_genome_reads/aligned_reads2"
    reference_genome="/projects/p31922/05_genome/polished_genomes/asm.np2.p_ctg.fa"

    # Loop through all forward FASTQ files in the input directory
    for forward in "$input_dir"/*.1.fq.gz; do
    # Extract the sample name (without path and forward extension)
    sample_name=$(basename "$forward" .1.fq.gz)

    # Define the corresponding reverse read file
    reverse="$input_dir/${sample_name}.2.fq.gz"

    # Check if the reverse file exists
    if [[ -f "$reverse" ]]; then
        # Run bwa mem for paired-end reads
        bwa mem "$reference_genome" "$forward" "$reverse" > "$output_dir/${sample_name}.sam"

        echo "Processed $sample_name"
    else
        echo "Warning: Reverse read file for $sample_name not found!"
    fi
    done    
    

Check quality of mapping on 1 or 2 individuals with code:

    samtools flagstat <SAMPLE.sam>
    
    
    Mapped a random indiv (WC8045):
    
    samtools flagstat WC8045.sam
    
    181791 + 0 in total (QC-passed reads + QC-failed reads)
    993 + 0 secondary
    0 + 0 supplementary
    0 + 0 duplicates
    144144 + 0 mapped (79.29% : N/A)
    180798 + 0 paired in sequencing
    90399 + 0 read1
    90399 + 0 read2
    129658 + 0 properly paired (71.71% : N/A)
    135958 + 0 with itself and mate mapped
    7193 + 0 singletons (3.98% : N/A)
    6176 + 0 with mate mapped to a different chr
    833 + 0 with mate mapped to a different chr (mapQ>=5)

On reads truncated to 100bp:
K=19 (default): 79.29% mapped
K=16: 75.66% mapped
K=14: 75.66
K=12: 76.03
K=10: 75.66
All are not very good, >90% is ideal - poorly sequenced indiv most likely - missing data issues again. We'll see if mapped alignment is any better. Explore removing individuals like this - very few starting reads (181k)

On untruncated reads:
K=19 (default): 77.90% mapped
K=16: 77.9% mapped
K=14: 77.9
Sticking with truncated reads. 


On individual PZ0010:
        
    (base) [amd9539@quser31 06_genome_reads]$ samtools flagstat PZ0010.K19.sam
    2932394 + 0 in total (QC-passed reads + QC-failed reads)
    20640 + 0 secondary
    0 + 0 supplementary
    0 + 0 duplicates
    2655355 + 0 mapped (90.55% : N/A)
    2911754 + 0 paired in sequencing
    1455877 + 0 read1
    1455877 + 0 read2
    0 + 0 properly paired (0.00% : N/A)
    2558288 + 0 with itself and mate mapped
    76427 + 0 singletons (2.62% : N/A)
    2295204 + 0 with mate mapped to a different chr
    20744 + 0 with mate mapped to a different chr (mapQ>=5)


## Assess read alignment for all reads

    ## With RADseq reads mapped to a reference genome, convert sam to bam, index, sort and get quality stats

    module load samtools

    # Specify the directory containing the .bam files
    DIRECTORY="/projects/p31922/06_genome_reads/aligned_reads/"

    # Navigate to the directory
    cd "$DIRECTORY" || { echo "Directory not found"; exit 1; }

    # Check if there are any .sam files  if ls *.sam 1> /dev/null 2>&1; then
    for file in *.sam; do
        echo "Processing $file..."

        # Convert SAM to BAM
        samtools view -bS "$file" -o "${file%.sam}.bam"

        # Check if the BAM conversion was successful
        if [[ -f "${file%.sam}.bam" ]]; then
            echo "Created ${file%.sam}.bam"

            # Sort the BAM file
            samtools sort "${file%.sam}.bam" -o "${file%.sam}_sorted.bam"
            echo "Sorted ${file%.sam}.bam to ${file%.sam}_sorted.bam"

            # Index the sorted BAM file
            samtools index "${file%.sam}_sorted.bam"
            echo "Indexed ${file%.sam}_sorted.bam"

            # Generate flag statistics
            samtools flagstat "${file%.sam}_sorted.bam" > "${file%.sam}_flagstat.txt"
            echo "Generated flag statistics for ${file%.sam}_sorted.bam"

        else
            echo "Failed to create ${file%.sam}.bam"
        fi
    done
    else
    echo "No .sam files found in the directory."
    fi


## Create summary flagstat to assess overall mappping quality

    # Set the directory containing the .flagstat.txt files
    FLAGSTAT_DIR="/projects/p31922/06_genome_reads/aligned_reads/"

    # Create a summary file
    SUMMARY_FILE="summary_flagstat.txt"
    echo -e "Sample\tQC_Passed_Reads\tMapped_Reads\tMapped_%\tProperly_Paired_Reads\tProperly_Paired_%" > "$SUMMARY_FILE"
    # Loop through all .flagstat.txt files
    for file in "$FLAGSTAT_DIR"/*_flagstat.txt; do
      if [[ -e "$file" ]]; then
    # Extract the sample name (remove directory and extension)
    SAMPLE=$(basename "$file" _flagstat.txt)

    # Extract relevant information
    QC_PASSED=$(grep "in total" "$file" | awk '{print $1}')
    MAPPED=$(grep "mapped (" "$file" | awk '{print $1}')
    MAPPED_PCT=$(grep "mapped (" "$file" | awk -F '[()]' '{print $2}')
    PROPERLY_PAIRED=$(grep "properly paired (" "$file" | awk '{print $1}')
    PROPERLY_PAIRED_PCT=$(grep "properly paired (" "$file" | awk -F '[()]' '{print $2}')

    # Append the data to the summary file
    echo -e "$SAMPLE\t$QC_PASSED\t$MAPPED\t$MAPPED_PCT\t$PROPERLY_PAIRED\t$PROPERLY_PAIRED_PCT" >> "$SUMMARY_FILE"
      else
    echo "No _flagstat.txt files found in $FLAGSTAT_DIR."
    exit 1
      fi
    done

    echo "Summary of all _flagstat.txt files saved to $SUMMARY_FILE."

    
### Export summary_file to R Studio...





# VCF tools filtering steps/notes
* Beginning with PO (population code) and following speciation genomics walkthrough to assess data - I moved all population.snps.vcf files out of their individual pop files into a single vcf folder and renamed as "PO.pop.snps.vcf"
*  https://speciationgenomics.github.io/filtering_vcfs/
*  https://github.com/cheyennelmoore/Erigenia-Analyses 
    * Move the folder with vcf files to home directory
        * `cp -r vcf /home/amd9539/vcf`
    * Create variables for easy coding (change PO to appropriate name)
        * `SUBSET_VCF=~/vcf/PO.pop.snps.vcf`
        * `OUT=~/vcf/PO_out/PO`
    * load vcf tools
        * `module load vcftools/0.1.17`
    * Calculate:
        * Allele frequency
            * `vcftools --gzvcf $SUBSET_VCF --freq2 --out $OUT --max-alleles 2`
        * Mean depth per indiv
            * `vcftools --gzvcf $SUBSET_VCF --depth --out $OUT`
        * Mean depth per site
            * `vcftools --gzvcf $SUBSET_VCF --site-mean-depth --out $OUT`
        * Site quality
            * `vcftools --gzvcf $SUBSET_VCF --site-quality --out $OUT`
        * Proportion missing data/indiv
            * `vcftools --gzvcf $SUBSET_VCF --missing-indv --out $OUT`
        * Proportion missing data/site
            * `vcftools --gzvcf $SUBSET_VCF --missing-site --out $OUT`
        * Het and inbred per indiv
            * `vcftools --gzvcf $SUBSET_VCF --het --out $OUT`
    * Move all outputs onto laptop and process in R Studio per tutorial to select filtering metrics
    * I have chosen and set as below, but first, set in and output file for filtered .vcf
        * `VCF_IN=~/vcf/FR.pop.snps.vcf`
        * `VCF_OUT=~/vcf/FR.pop.snps.filtered.vcf`
    * Filtering metrics:
        * `MAF=0.2`
        * `MISS=0.75`
        * `MIN_DEPTH=10`
        * `MAX_DEPTH=50`
    * Filter with:
        * vcftools --gzvcf $VCF_IN --remove-indv LO0003 --remove-indels --maf $MAF --max-missing $MISS --min-meanDP $MIN_DEPTH --max-meanDP $MAX_DEPTH --minDP $MIN_DEPTH --maxDP $MAX_DEPTH --recode --stdout | gzip -c > $VCF_OUT

  



# Workflow for VCF Tools and Plink to set up PCA in R, view missing data, and view population stats

* Go to parameter folder that contains assembly data and populations.snps.vcf file
    * create plink folder
        * mkdir plink 
    * add scaffold_ to snp ID in new vcf file within plink folder so plink doesn't think it's human chromosome data and give you errors (https://groups.google.com/g/plink2-users/c/b5bP9YZfkFU)
        * `sed -E 's/^([0-9])/scaffold_\1/g' populations.snps.vcf > plink/pop.snps2.vcf `
    * go into plink folder and load VCF tools
        * `cd plink`
        * `module load vcftools/0.1.17`
    * create plink file from vcf
        * `vcftools --vcf pop.snps2.vcf --plink --out plink` 
    * load plink (https://speciationgenomics.github.io/pca/ <- great guide from here on with more details)
        * `module load plink/1.9`
    * perform linkage pruning - i.e. identify prune sites 
        * `plink --vcf pop.snps2.vcf --double-id --allow-extra-chr --set-missing-var-ids @:# --indep-pairwise 50 10 0.1 --out kplink`
    * create pca eigenvalue files
        * `plink --vcf pop.snps2.vcf --double-id --allow-extra-chr --set-missing-var-ids @:# --extract kplink.prune.in --make-bed --pca --out kplink$$$`
    * Export the .eigenval and .eigenvec files to your laptop and use to create a PCA in R Studio...
    * 
* In R Studio... 
    * '

* Check missing data
    * `vcftools --vcf pop.snps2.vcf --missing-indv --out MissingData_322$$$.imiss`
    * Create out.imiss file
        * Move to laptop and view in excel to sort by missing data

# 4/24/23 Using VCF tools and Plink to filter data and cluster in a PCA

* For help:
    * VCFtools: https://vcftools.github.io/man_latest.html
    * Plink: https://www.cog-genomics.org/plink/2.0/general_usage
* Both vcf tools and Plink are installed on Quest, load software with: 
    * module load vcftools/0.1.17
    * module load plink/1.9
        * plink/2.0 is listed as preferred but according to the plink manual, plink2 is still undergoing alpha and beta testing, so plink1.9 should be utilized primarily - plink2 should be used as a 'supplement' for specific tasks. 
* STEP 1 - create plink file from VCF file if not done during denovo assembly (need to create plink folder (mkdir) within parameter folder you are using):
    * vcftools --vcf [input folder/populations.snps.vcf] --plink --out [output folder/plink]
        * this took 2 seconds and 1.5MB of memory
* STEP 2 - TROUBLESHOOTING
    * Error: Invalid chromosome code 27 resolution discussion - https://groups.google.com/g/plink2-users/c/b5bP9YZfkFU
    * Terminology - contig is a continuous stretch of nucleotide sequences while scaffold is a portion of the genome consisting of contigs and gaps. Both contig and scaffold are reconstructed genomic sequences
    

# 4/17/23 Compiling stats for each denovo assembly parameter

* In each parameter folder (322, 333 etc.) open the populations.log file to find loci kept and loci variants
* Enter into excel sheet to compare and make decisions based on loci 
* Should I assemble with no population map to generate stats for the subsample?
    * Yes - run wit and without pop map and we will run PCA
* Use plink to generate PCA for assembly run on each population with each sample 
* THEN use PCA to assess if sublibrary or pop better explain clustering

# Conda location


  environment location: /home/amd9539/.conda/envs/STACKS
  
# 3/15/23 running STACKS with Slurm via notepad ++

## Create script in Notepadd ++ that contains code for STACKS programing

QUESTIONS
* How do I view depths of coverage output? What is the filename? A: gstacks.log.distribs
* Does specifying the number of cores in stacks denovo assembler change how slurm allocates cpu usage in quest?
* Submitting jobs successfully to slurm but commands not actually running - is STACKS not uploaded properly/not being called in properly? No errors in slurm error logs, nothing in slurm output folder.

NOTES
* To save file as ''.sh', Save As, select 'All Types', and end file name with .sh
* NEED TO RUN STACKS WITH parameters: m,M,n; 3,2,2; 4,2,2; 5,2,2; 3,3,3; 3,4,4; 3,5,5
* Use STACKS on QUEST - called with module load stacks/2.64-gcc-9.2.0
* JOB 6421078 COMPLETED (3/27) - "denovo_map.pl" on 49 samples, QUEST allocation used: 1 node, 4 cores/node, Time = 2hr:30min, Memory = 3GB (does memory increase with more data? NO, but add a couple of GB f it's close to request)
* Scaling up (3/28) - 60 samples randomly selected from all pools/sublibraries (23 from 1 & 2, 14 from 3; sublibrary 3 has only 30 samples). Script set to run on 6 different parameter sets = 2hr30min x 6 = 15hr, +10 samples => use 24 hr to be safe?? Memory to use?
    * submitted job with 24hr time, 4GB memory, 1 node, name = Test_Stacks (remember to change job name next time)
        * Job cancelled - 
            ID: 6536925
            Cluster: quest
            User/Group: amd9539/amd9539
            State: NODE_FAIL (exit code 1) => ????
            Nodes: 1
            Cores per node: 4
            CPU Utilized: 00:00:00
            CPU Efficiency: 0.00% of 00:12:24 core-walltime
            Job Wall-clock time: 00:03:06
            Memory Utilized: 0.00 MB (estimated maximum)
            Memory Efficiency: 0.00% of 6.00 GB (6.00 GB/node)


## Create script in Notepad ++ for SLURM that references STACKS script as .sh file with bash




  
  