Several script useful to analyse or explore genomic data
git clone https://github.com/Yuma248/MELFUtilities
This scritp extracts sequences from a dDocent reference fasta file, according to infromation from a VCF file with your filtered SNPs.
Usage:
ExtSeq.pl
-inR <path to reference file fron dDocent, with contig sequences>
-inV <path to VCF file with SNP>
-out <output fasta file name>
Example:
ExtSeq.pl -inR /dDocent/reference.fasta -inV /dDocent/Filtered_SNPs.vcf -out myseq.fasta
This script will extract the best quality SNP, first SNP, less missing data SNP and a random SNP per contig. It is usefull when SNPs were call de novo to reduce linkage, by selecting one SNP per contig.
Usage:
ExtractBQFRLMRM.pl
-inf <input vcf_file>
-bq <output best quality SNP, default Best.vcf>
-frt <First SNP ouput, default First.vcf>
-rndm <random SNP output, default Random.vcf>
-lm <output SNP with lest missing data, defaul Lmissin.vcf>
For example:
ExtractBQFRLMRM.pl -inf /my/input.vcf -bq /my/bestSNP.vcf -frt /my/firstSNP.vcf -rndm /my/randomSNP.vcf -lm /my/lessmissingSNP.vcf
This script will upload a file to a figshare repository, it needs the whole path of the file to be uploaded, a URL from figshare, and a token to access the repository or project. If the file is not required to be uploaded in a specific project or item, you can use the general URL https://api.figshare.com/v2/account/articles. But if you need to put in a specific project or item (folder) you will require the project ID (using the -p option) and/or item ID (usign the -f option). If you already upload one file in the project and item, you can integrate the item ID in the URL and add just /files at the end, see the example.
Usage:
UpFigshare.sh
-i <path to input file>
-u <URL of figshare account/project/item, required if you have specific project and item ready, default https://api.figshare.com/v2/account/articles\>
-t <Access token for the figshare account>
-p <project ID, this is required just if you require to upload the file in specific project>
-f <item ID, this is required just if you need to upload the file in specific item (folder)>
-n <item name, this is required if you are creating a new item (folder) and you want a specific name, default NEW UPLOAD>
For example:
UpFigshare.sh -i /sanpper/genome/genomev2.fasta -u https://api.figshare.com/v2/account/articles/ -t 75050303931z87ab7c72038ab9eaf02d853766bd8f7cc695f390d5b9cdeda1fd230c462a7cb7c7a67ef7c507f27f3fde647c4145667664533374d54609ef477874c4aa11 -p 24892
This script performs basic filtering on a dart CSV file and output a table with the remaining SNPs after each filtering step, the basic stats (#samples, #loc. #polymorphic loc, observerd heterozygozity, espected heterozygozity, and Fis inbreeding coefficient) after filtering, a FST plot, a PCA plot and a Admixture/Structure plot.
To use this fuction in R you will need the to install and load the next R libraries:
library(dartR)
library(tidyverse)
library(vegan)
library(related)
library(igraph)
library(reshape2)
library(ggplot2)
library(HardyWeinberg)
library(ggpubr)
library(gridExtra)
library(LEA)
And either copy and save the fuction in R or use the command
source("path/to/MELFUtilties/basic_filtering_dart.R")
You will have to upload your dart data and samples metadata (sample pop map or pop info file) in R, using commands like:
dartData <- gl.read.dart(filename="Report_Sps0000_SNP_mapping.csv",ind.metafile="SamplesAllMetadata.csv")
meta <- read.csv("SamplesAllMetadata.csv", header=TRUE)
Usage:
basic_filter(dartData, "output prefix", maxmisi = 50, mincalL = 0.80, minrep = 0.99, minmaf = 0.03, secd = TRUE, HWEF = TRUE, depthr = c(5,75), npopsHE = 5, maxsim = 0.85, PDFplots = TRUE)
dartD - DArT genotype file
name - prefix for output files (string)
maxmisi - maximum proportion of missing data per individual
mincalL - minimum call rate per locus, inverse of maximum missing data por locus
minrep - minimum reproducibility
minmaf - minimum minor allele frequency
maxsim - max reletedness to remove duplicates
secd - remove secundaries (TRUE, FALSE)
depthr - min and maximum depth of coverage
HW - filter for HWE (TRUE, FALSE)
npopsHE - Maximum number of population where a locus can be out of HWE
PDFplots - If you want to create a PDF with basic results (filter steps table, diversity table, Fst heat map, PCA)
Example:
filterdata <- basic_filter(YumadartData, "Yuma", maxmisi = 50, mincalL = 0.80, minrep = 0.99, minmaf = 0.03, secd = TRUE, HWEF = TRUE, depthr = c(5,75), npopsHE = 5, maxsim = 0.85, PDFplots = TRUE)