ALGINMENT2VCF

First Download Brachypodium Genome reference file

https://www.ncbi.nlm.nih.gov/guide/howto/dwn-genome/

Select Genomes FTP site .....

Select /refseq index

Select /plant Index

Select /Brachypodium _distachyon

Select /all_assembly_versions

I choose version one v1

I choose an ascension run ERR4835478 and used fastq-dump ERR4835478 to generate a fastq file

fastq-dump ERR4835478

I use BWA to align Fasta with FASTQ and generate BAM file

~9 Hours @ 18.1 GB

Take a look at the HEader of The FASTQ file

We run fastqc R package on our FASTQ file with the following R script

setwd("/home/michael/Desktop/Alignment")
install.packages("fastqcr")
library(fastqcr)
fastqc_install()
fastqc(fq.dir = "/home/michael/Desktop/Alignment/FASTQ",threads = 4)
qc.file <- "/home/michael/Desktop/Alignment/FASTQ/FASTQC/ERR4835478_fastqc.zip"
qc <- qc_read(qc.file)
names(qc)
#These commands are only necessary to plot in RStudio Viewer, otherwise the most important plots are automatically genearted to HTML file in FASTQC directory

qc_plot(qc, "summary")
qc_plot(qc, "Basic statistics")
qc_plot(qc, "Per base sequence quality")
qc_plot(qc, "Per base sequence content")
qc_plot(qc, "Per sequence GC content")
qc_plot(qc, "Per sequence quality scores")
qc_plot(qc, "Sequence duplication levels")
qc_plot(qc, "Per base N content")
qc_plot(qc, "Sequence length distribution")
qc_plot(qc, "Sequence duplication levels")
qc_plot(qc, "Overrepresented sequences")
qc_plot(qc, "Adapter content")
qc_plot(qc, "Kmer content")

Open the HTML file in the FASTQC directory Generated by fastqc function

I have taken several screenshots

Check the Data Structure

Index your reference genome FASTA file first

Use samtools to sort BAM file

Index the Sorted Reads

Make a VCF file with BAM, Indexed Bam, and Reference Genome FASTA

At this point you should have 6 directories and 17 files

Rename the sample to match Ascension Identifier

VCF Header and Info Fields

Use ChromQual to plot important figures with R Script

setwd("/home/michael/Desktop/Alignment/plots")
library(QTLseqr)
library(tinytex)
library(vcfR)
library(tidyr)
library(ggplot2)
library(dplyr)
library(ggrepel)
library(ggpubr)
library(data.table)
# Define a vector of Chromosomes
Chroms <- c("NC_016131.3","NC_016132.3","NC_016133.3","NC_016134.3","NC_016135.3")
file <- "vcfnewsamplename.vcf.gz"
QTLseqr::ChromQual(file = file,chromlist = Chroms, windowSize = 1e+05, HighLimQuality = 15000, scalar = 1, ncol = 5, binwidth1 = 10, binwidth2 = 100, p1 = TRUE, p2 = TRUE, p3 = TRUE, p4 = TRUE, p5 = TRUE)

Standard R Console Output

And Plots p1, p2, and p5 as Generated in View Panel

Using RMVP Package plot SNP densities per chromosome/contig

library(rMVP)
sample<-"ERR4835478"
pathtosample <- "/home/michael/Desktop/Alignment/plots/vcfnewsamplename.vcf"
out<- paste0("mvp.",sample,".vcf")
memo<-paste0(sample)
dffile<-paste0("mvp.",sample,".vcf.geno.map")

message("Making MVP data S1")
MVP.Data(fileVCF=pathtosample,
         #filePhe="Phenotype.txt",
         fileKin=FALSE,
         filePC=FALSE,
         out=out
)
message("Reading MVP Data S1")
df <- read.table(file = dffile, header=TRUE)
message("Making SNP Density Plots")
MVP.Report.Density(df[,c(1:3)], bin.size = 100000, col = c("blue", "yellow", "red"), memo = memo, file.type = "jpg", dpi=300)