Building a bioinformatic analysis platform on Windows 10
The Windows Subsystem for Linux (WSL) allows users to install Linux subsystems directly on a Windows 10 system. It can easily run Linux commands and install Linux software to avoid the installation of third-party virtual machine software. The advantage of WSL is that it makes better use of computer memory and does not require copying files between the host and the virtual machine.
Configuration of WSL
Timing ~1 min
System requirements: Windows 10 Version 1709, Build 16299, or above 64-bit systems.
1. Enable WSL: Open “Settings”, click “Apps”, then find and click “Programs and Features”, click “Turn Windows features on or off”, find “Windows Subsystem for Linux” and check the box, click “OK”, and restart the computer (Supplementary Video 1).
Timing ~59 min
2. Open the Microsoft Store, search Ubuntu, and choose to install Ubuntu 18.04 LTS. Follow the prompts to set up your username and password. Here, we create an account with the username “bio” (Supplementary Video 2). When the installation is finished, we need to do some configuration on the system (Supplementary Video 3).
3. Enter the following command in the terminal to update the source:
$sudo apt-get update
4. Set the password for root.
$sudo passwd root
5. Enable the CUDA-aware MPI.
For Linux 64, the CUDA awareness support may be disabled by default. Users should enable the support by setting the environment variable to use OpenMPI. Check whether CUDA awareness support is enabled in the environment variable configuration file (~/.bashrc). If it is not enabled, enter the following commands in the terminal.
$echo OMPI_MCA_opal_cuda_support=true >> ~/.bashrc
6. Installation of Miniconda
(A) Installation of Miniconda on Linux
(i) Here, Miniconda will be installed; go to the official website, and select the installation file suitable for your system and Python version (Supplementary Video 4).
(ii) Start installation
Keep pressing “Enter” key when prompted to visualize the license agreement, enter “yes” and press “Enter” to continue. Press “Enter” to confirm the default installation location. Miniconda is installed in the miniconda3 directory, under the user’s home directory. Type “yes” and press “Enter” to initialize miniconda3. Finally, type the command “source ~/.bashrc” in the terminal.
(iii) Set up the Bioconda channel. Add the channels by entering the following three commands in the terminal.
$conda config --add channels defaults
$conda config --add channels bioconda
$conda config --add channels conda-forge
(B) Install Miniconda on MacOS
(i) Installation of Miniconda3
(ii) Add channels of Bioconda
$conda config --add channels defaults
$conda config --add channels bioconda
$conda config --add channels conda-forge
Installation of PGCGAP (Supplementary Video 5).
Timing ~34 min
7. Create a pgcgap environment and install PGCGAP (Usually specify the latest version of PGCGAP).
$conda install -y mamba
$mamba create -n pgcgap pgcgap=1.0.34 java-jdk=8.0.112
8. Activate the pgcgap environment.
$conda activate pgcgap
9. Set up the COG database.
10. Exit the pgcgap environment.
Step by Step examples
Timing ~20 h
The usage and parameters of PGCGAP can be viewed by typing “pgcgap -h” in the terminal. Next, we show how to run all the modules of the PGCGAP through a dataset.
11. Download and decompress the example dataset.
$tar -zxvf PGCGAP_Examples.tar.gz
In this example, the working directory is located at the H drive. All hard disks in Windows were mounted in the “/mnt” directory of Ubuntu Linux. The “PGCGAP_Examples/Reads/Illumia” directory contains six Illumina Hiseq paired-end reads of Escherichia coli; the “PGCGAP_Examples/Reads/Oxford” directory contains the Oxford Nanopore reads of Escherichia coli K12; and the “PGCGAP_Examples/Reads/PacBio” directory contains the Pacific Biosciences released P6-C4 chemistry reads of Escherichia coli K12. “PGCGAP_Examples/Reads/MG1655.gbff” is the GenBank format file of E. coli K-12 substr. MG1655, and will be used as the reference genome. The “PGCGAP_Examples/Reads/Hybrid” directory contains two short reads files and one long reads file of the same strain. “PGCGAP_Examples/Other_inputs/ proteins.fas” contains 18 protein sequences of MFS transporter from several bacterial species.
12. Activate the pgcgap environment.
$conda activate pgcgap
13. Example 1: Genome assembly with Illumina reads.
Paired-end reads of six strains in the directory “Reads/Illumina/” are used as inputs. In the dataset, the naming format of the genome is “strain_1.fastq.gz, ” and “strain_2.fastq.gz”. The string after the strain name is “_1.fastq.gz”, and its length is 11, so “--suffix_len” was set to 11. Users can choose “abyss”, “spades”, and “auto” for genome assembly. The assembly speed with “abyss” is faster, and the assembly quality with “spades” is better. Taking into account the speed and quality of assembly, we suggest using the “auto” mode for assembly. “--filter_length” is set here to remove sequences shorter than 200 bp from the assembled genomes.
$pgcgap --Assemble --platform illumina --assembler abyss --filter_length 200 --ReadsPath Reads/Illumina --reads1 _1.fastq.gz --reads2 _2.fastq.gz --kmmer 81 --threads 4 --suffix_len 11
$pgcgap --Assemble --platform illumina --assembler spades --filter_length 200 --ReadsPath Reads/Illumina --reads1 _1.fastq.gz --reads2 _2.fastq.gz --threads 4 --suffix_len 11
$pgcgap --Assemble --platform illumina --assembler auto --filter_length 200 --ReadsPath Reads/Illumina --reads1 _1.fastq.gz --reads2 _2.fastq.gz --kmmer 81 --threads 4 --suffix_len 11
New directories and documents are generated after the program is completed. The assembly results for each genome are in the “Results/Assembles/Illumina” directory, while all scaffolds of the strains are stored in “Results/Assembles/Scaf/Illumina”. “*.filtered.fas” is the genome with short sequences removed. “*.prefilter.stats” describes the status of the genome before filtering, and “*.filtered.stats” describes the status of the genome after short sequence filtering. While “abyss” was chosen as the assembler, users are advised to check the assembly stats file (such as Results/Assembles/Illumina/SRR9620252_assembly/SRR9620252-stats.tab) of each genome to ensure that the value of N50 is greater than 50,000 bp. The file “scaf.list” under the working directory contains the absolute path of all genomes, which will be used as the input file of module ANI.
14. Example 2: ONT reads assembly and polishing.
The Oxford nanopore produces only one read file (“Reads/Oxford/oxford.fasta”), so only the parameter of “--reads1” needs to be set. Here, the value “.fasta”. “--genomeSize” is the estimated genome size, and users can check the genome size of similar strains in the NCBI database for reference. The parameter is set to “4.8m”. The suffix of the reads file here is “.fasta” and its length is 6, so “--suffix_len” is set to 6.
$pgcgap --Assemble --platform oxford --ReadsPath Reads/Oxford --reads1 .fasta --genomeSize 4.8m --threads 4 --suffix_len 6 --filter_length 200
The results are stored in the “Results/Assembles/Oxford” and “Results/Assembles/Scaf/Oxford” directories. The former contains all intermediate files and genome files, while the latter contains only the assembled genome.
15. Example 3: PacBio reads assembly and polishing.
PacBio also produces only one read file (“Reads/PacBio/pacbio.fastq”); the parameter settings are similar to those of Oxford. The strain name is “pacbio” with the suffix “.fastq” and the suffix length is 6, so “--suffix_len” was set to 6.
$pgcgap --Assemble --platform pacbio --ReadsPath Reads/PacBio --reads1 .fastq --genomeSize 4.8m --threads 4 --suffix_len 6 --filter_length 200
The results are stored in the “Results/Assembles/PacBio” and “Results/Assembles/Scaf/PacBio” directories. The former contains all intermediate files and genome files, while the latter contains only the assembled genome.
16. Example 4: hybrid assembly of short reads and long reads and polishing.
Paired-end short reads and long reads in the directory “Reads/Hybrid/” are used as inputs. Illumina reads and long reads had been obtained from the same isolates.
$pgcgap --Assemble --platform hybrid --ReadsPath Reads/Hybrid --short1 short_reads_1.fastq.gz --short2 short_reads_2.fastq.gz --long long_reads_high_depth.fastq.gz --threads 4
The results are stored in the “Results/Assembles/Hybrid” directory, and the final assembly is named “assembly.fasta”.
17. Example 5: Gene prediction and annotation.
Here, the assembly results of Illumina reads are taken as inputs (“Results/Assembles/Scaf/Illumina/*.filtered.fas”). The suffix of the genome is “.filtered.fas”. When running the program, the value of the “--Scaf_suffix” parameter cannot be quoted. Here, .filtered.fas should not be quoted.
$pgcgap --Annotate --scafPath Results/Assembles/Scaf/Illumina --Scaf_suffix .filtered.fas --genus Escherichia --species coli --codon 11 --threads 4
The generated files are stored in the “Results/Annotations” directory, and files in the directories “Results/Annotations/AAs”, “Results/Annotations/CDs” and “Results/Annotations/GFF” will be used for subsequent analysis.
18. Example 6: Constructing the single-copy core protein tree and core SNP tree.
The phylogenetic trees of single-copy core proteins and single-copy core gene SNPs will be constructed using the six E. coli genomes sequenced by Illumina as datasets. The input files are the amino acid sequence files (“Results/Annotations/AAs/*.faa”) and the nucleotide sequence files (“Results/Annotations/CDs/*.ffn”) obtained by genome annotation. Amino acid files and nucleotide files must be suffixed with “.faa” and “.ffn”, respectively. The “.faa” and “.ffn” files of the same strain should have the same prefix name. The name of protein IDs and gene IDs in the amino acid file and nucleotide file should start with the strain name. The Prokka12 software was suggested to generate the input files. “strain_num” is a very important parameter to specify the number of genomes for analysis. Three parameters including --fasttree, --bsnum (default), and --fastboot are provided to choose the method for phylogenetic tree constructing.
$pgcgap --CoreTree --CDsPath Results/Annotations/CDs --AAsPath Results/Annotations/AAs --codon 11 --strain_num 6 --threads 4
The result files are stored in the directory “Results/CoreTrees”. “ALL.core.protein.fasta.gb.treefile” and “ALL.core.snp.fasta.gb.treefile” are the phylogenetic tree files constructed based on the aligned sequences of the single-copy core proteins and the core SNPs with the best-fit model of evolution, respectively. The best-fit model of evolution for the protein alignments and SNPs alignments can be found in the file “ALL.core.protein.fasta.gb.iqtree” and “ALL.core.snp.fasta.gb.iqtree”, respectively. While --fasttree was sprcified, “ALL.core.protein.nwk” will be outputted as the phylogenetic tree of single-copy core proteins. Users can import tree files into MEGA30 or iTOL31 to view the topology.
19. Example 7: Constructing the single-copy core protein tree only.
If the “--CDsPath” was set to “NO”, the nucleotide files will not be needed, and the phylogenetic tree of core SNPs will not be constructed.
$pgcgap --CoreTree --CDsPath NO --AAsPath Results/Annotations/AAs --codon 11 --strain_num 6 --threads 4
20. Example 8: pan-genome analysis and phylogenetic tree construction.
(A) GFF3, also known as generic feature format version 3, is a tab-delimited, plain text file used to describe features of DNA, RNA, and protein sequences. For Pan-genome analysis, GFF3 files (With “.gff” as the suffix) of each strain are placed into a directory (“Results/Annotations/GFF/*.gff”). We strongly recommend using Prokka12 to generate the aforementioned files. If the Annotate module runs first, the files are automatically generated. Users can set the minimum percentage identity for blastp by the parameter “identi”, and only an integer number is allowed. Users can choose any of the three parameters including --fasttree, --bsnum (default), and --fastboot to construct the phylogenetic tree.
$pgcgap --Pan --codon 11 --strain_num 6 --threads 4 --GffPath Results/Annotations/GFF --PanTree --AAsPath Results/Annotations/AAs
The results are stored in the “Results/PanGenome” directory. A spreadsheet named “gene_presence_absence.csv” lists each gene and which samples contained it. Users can take the gene_presence_absence.csv file and a trait file to conduct pan-genome wide association studies with the Scoary34 software. At the same time, some visual results (“*.pdf”) are also outputted. A phylogenetic tree based on the aligned sequences of single-copy core proteins called by Roary will be constructed automatically if the parameter “PanTree” was provided. The tree file can be found in the directory “Results/PanGenome/Core/” with a name of “Roary.core.protein.fasta.gb.treefile” or “Roary.core.protein.fasta.gb.nwk”. The best-fit model of evolution for the protein alignments can be found in the file “Roary.core.protein.fasta.gb.iqtree”.
21. Example 9: Inference of orthologous gene groups.
The input files are also the amino acid sequence files suffixed with “.faa” (“Results/Annotations/AAs/*.faa”). Users can choose any of the three parameters including --fasttree, --bsnum (default), and --fastboot to construct the phylogenetic tree.
$pgcgap --OrthoF --threads 4 --AAsPath Results/Annotations/AAs
The resulting files are placed in the “Results/OrthoFinder/Results_orthoF” directory. The file “Single.Copy.Orthologue.fasta.gb.treefile” in directory “Single_Copy_Orthologue_Tree” is the tree file of single-copy orthologues, and the best fit module can be found in the file “Single.Copy.Orthologue.fasta.gb.iqtree” under the same directory.
22. Example 10: Compute whole-genome Average Nucleotide Identity.
The input file named “scaf.list” contains the absolute path of each genome, one per line. If the “--Assemble” function is run first, the list file is generated automatically. The value of the parameter “--Scaf_suffix” depends on the actual situation, here is “.filtered.fas”.
$pgcgap --ANI --threads 4 --queryL scaf.list --refL scaf.list --ANIO Results/ANI/ANIs --Scaf_suffix .filtered.fas
The results are stored in the “Results/ANI” directory. The file “ANI” contains comparison information of genome pairs. The document is composed of five columns, each of which represents the query genome, reference genome, ANI value, count of bidirectional fragment mappings, and total query fragments. A heat map file “ANI_matrix.pdf” is generated.
23. Example 11: Genome and metagenome similarity estimation using MinHash
This requires genome files (complete or draft) in a directory as inputs (Default: Results/Assembles/Scaf/Illumina).
$pgcgap --MASH --scafPath Results/Assembles/Scaf/Illumina --Scaf_suffix .filtered.fas
The results are stored in the “Results/MASH” directory. The file “MASH” shows the pairwise distance between pair genomes, and each column represents Reference-ID, Query-ID, Mash-distance, P-value, and Matching-hashes. A heat map file named “MASH_matrix.pdf” is generated to describe the similarity of each genome pair.
24. Example 12: COG annotation.
The input files are also the amino acid sequence files suffixed with “.faa” (“Results/Annotations/AAs/*.faa”).
$pgcgap --pCOG --threads 4 --strain_num 6 --AAsPath Results/Annotations/AAs
The results are stored in the “Results/COG” directory. The super COG table of each strain (“*.2Scog.table”) and its plot (“*.2Scog.table.pdf”) will be generated. “All_flags_relative_abundances.table” is a table containing the relative abundance of each flag for all strains, while “All_flags_relative_abundances.pdf” is the corresponding visualization result.
25. Example 13: Variants calling and phylogenetic tree construction based on a reference genome.
The six genomes sequenced by Illumina were chosen as datasets (“Reads/Illumina/*.gz”). Escherichia coli K-12 substr. MG1655 was selected as the reference genome and the reference file “MG1655.gbff” in the GenBank format is stored in the “Reads” directory. The absolute path of the reference genome (here is “/mnt/h/PGCGAP_Examples/Reads/MG1655.gbff”) is required to run the program.
$pgcgap --VAR --threads 4 --refgbk /mnt/h/PGCGAP_Examples/Reads/MG1655.gbff --ReadsPath Reads/Illumina --reads1 _1.fastq.gz --reads2 _2.fastq.gz --suffix_len 11 --strain_num 6 --qualtype sanger
The resulting files are stored in the “Results/Variants” directory, where the “Core” directory contains the core SNPs of all strains and their phylogenetic tree.
26. Example 14: Screening of contigs for antimicrobial and virulence genes
This requires genome files (complete or draft) in a directory as inputs (Default: Results/Assembles/Scaf/Illumina). Users can choose one of the following databases for analysis: argannot32, card33, ecoh34, ecoli_vf (https://github.com/phac-nml/ecoli_vf), ncbi35, plasmidfinder36, resfinder37, MEGARes38, and vfdb39. Users can set the value of parameter “db” to “all” to search against to all of the databases.
$pgcgap --AntiRes --scafPath Results/Assembles/Scaf/Illumina --Scaf_suffix .filtered.fas --threads 4 --db all --identity 75 --coverage 50
The resulting files are stored in the “Results/AntiRes” directory. “*.tab” files are screening results of each strain, and the “summary.txt” file contains a matrix of gene presence/absence for all strains.
27. Example 15: Perform all analyses for paired-end reads.
In fact, we provide many one-line commands to complete analyses of each module with original reads, which have been demonstrated in the examples of module Annotate, pCOG, CoreTree, Pan, OrthoF, ANI, MASH, and AntiRes. In addition, users can complete almost all analyses with original reads by just one line of commands through the module All. In short, the module All will call modules Assemble, Annotate, pCOG, CoreTree, Pan, OrthoF, ANI, MASH, and AntiRes one by one to perform the analyses. For the sake of flexibility, the module VAR was not integrated into module All and needed to be added when needed. Only the reads and reference files should be provided, and important parameters should also be provided.
$pgcgap --All --platform illumina --filter_length 200 --ReadsPath Reads/Illumina --reads1 _1.fastq.gz --reads2 _2.fastq.gz --suffix_len 11 --kmmer 81 --genus Escherichia --species coli --codon 11 --strain_num 6 --threads 4 --VAR --refgbk /mnt/h/PGCGAP_Examples/Reads/MG1655.gbff --qualtype sanger --PanTree
28. Example 16: Filter short sequences in the genome and assess the status of the genome.
“Assess” takes the assembled genomes as inputs. First, it assesses the stats of the genome; second, the sequences shorter than “--filter_length” are deleted from the genome. Finally, the stats of the filtered genome are assessed. The results files are stored in the same directory as the inputs.
$pgcgap --ACC --Assess --scafPath Results/Assembles/Scaf/Illumina --Scaf_suffix -8.fa --filter_length 200
29. Example 17: Construct a phylogenetic tree based on multiple FASTA sequences in one file.
“STREE” takes the file containing multiple-FASTA sequences as input. The parameter “--bsnum” represents the number of bootstraps.
$pgcgap --STREE --seqfile Other_inputs/proteins.fas --seqtype p --bsnum 500 --threads 4
The result files will be stored in the directory “Results/STREE”. The file “proteins.fas.aln.gb.treefile” is the final phylogenetic tree file. The best-fit model of evolution for the protein alignments can be found in the file “proteins.fas.aln.gb.iqtree”.