Mercurial > repos > rliterman > csp2

comparison CSP2/README.md @ 0:01431fa12065

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload"
author rliterman
date Mon, 02 Dec 2024 10:40:55 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:01431fa12065
1 <p align="center">
2 <img src="img/Temp_Logo.jpg" alt="drawing" width="400"/>
3 </p>
4
5 # CFSAN SNP Pipeline 2 (CSP2)
6 #### Dr. Robert Literman
7 **Office of Analytics and Outreach**
8 **Center for Food Safety and Applied Nutrition**
9 **US Food and Drug Administration**
10
11 **Current Release**: [v.0.9.5 (Oct-17-2024)](https://github.com/CFSAN-Biostatistics/CSP2/releases/tag/v.0.9.5)
12 **Last Push**: Oct-17-2024
13
14 **Important Note:** *CSP2 is currently under development, and has not been validated for non-research purposes. Current workflows and data processing parameters may change prior to full release version.*
15
16 ## CSP2 is a Nextflow pipeline for rapid, accurate SNP distance estimation from assembly data
17
18 CSP2 runs on Unix, with the handful of dependencies listed in the [Software Dependencies](#software-dependencies) section. CSP2 was developed to **(1)** improve on the speed of the [CFSAN SNP Pipeline (CSP)](https://peerj.com/articles/cs-20/?report=reader), **(2)** to reduce computational burden when analyzing larger isolate clusters, and **(3)** to remove the dependency for raw Illumina data. CSP2 relies on the accurate and rapid alignment of genome assemblies provided by [MUmmer](https://github.com/mummer4/mummer), which typically complete within seconds. This provides significant reductions in runtime compared to methods that rely on read mapping. The use of assemblies in place of sequencing data also means that:
19
20 - the amount storage needed can be substantially reduced,
21 - significantly less computational resources are required,
22 - as long as assemblies are available, isolates can be compared regardless of sequencing platform or whether publicly available sequence data even exists
23
24 CSP2 runs are managed via Nextflow, providing the user with an array of [customizations](#tips-for-configuring-csp2) while also facilitating module development and additions in future releases.
25
26 **Important Note**: *The software continues to be focused on the analysis of groups of bacterial genomes with limited evolutionary differences (<1000 SNPs). Testing is underway to determine how the underlying cluster diversity impacts distances estimates.*
27
28 ### CSP2 has two main run modes (See [Examples](#examples)):
29
30 #### 1) "Screening Mode" (*--runmode screen*): Used to determine whether query isolates are close to a set of reference isolates (e.g., lab control strains, strains related to an outbreak, etc.)
31 Given one or more user-provided reference isolates (*--ref_reads*; *--ref_fasta*), get alignment statistics and SNP distances between all reference and query isolates (*--reads*; *--fasta*)
32
33 #### 2) "SNP Pipeline Mode" (*--runmode snp*): Used to generate pairwise distances and alignments for a set of query isolates
34 Generate pairwise SNP distances and alignments for 2+ isolates (*--reads*; *--fasta*) based on comparisons to:
35 - One or more user-provided references (*--ref_reads*; *--ref_fasta*), or
36 - One or more reference isolates selected by RefChooser (*--n_ref*)
37
38 All CSP2 sequence comparisons happen at the assembly level, but if reads are provided CSP2 will perform a genome assembly using *SKESA*. In either case, CSP2 then calls MUMmer for alignment. If a sufficient portion of the reference genome is aligned (*--min_cov*), that data is passed through a set of filters that largely mimic those from the CFSAN SNP Pipeline, including the automated removal of:
39 - Sites from short alignments (*--min_len*)
40 - Sites from poorly aligned contigs (*--min_iden*)
41 - Sites close to the contig edge (*--query_edge*/*--ref_edge*)
42 - Sites from regions of high SNP density (*--dwin*/*--wsnps*)
43 - Multiply aligned sites
44 - Non-base sites (e.g., 'N' or '?')
45 - Heterozygous sites
46 - Indels (**for now**)
47
48 This final dataset is summarized into a *.snpdiffs* file, which contains:
49 1. A one-line header with alignment statistics
50 2. A BED file of contig mappings that pass QC
51 3. Information about SNPs (if present)
52
53 To avoid unnecessary realignment, once a .snpdiffs file is generated under a particular set of QC parameters (which is hardcoded into the .snpdiffs file as the "QC_String") these files can be used in other CSP2 runs via the *--snpdiffs* argument (if using the same QC parameters).
54
55 ---
56
57 ## Software Dependencies
58 The following software are required to run CSP2. Software version used during CSP2 development noted in parentheses.
59
60 - [Nextflow](https://www.nextflow.io/docs/latest/getstarted.html) (22.10.7)
61 - [Python](https://www.python.org/downloads/) (3.8.1)
62 - [pybedtools](https://pypi.org/project/pybedtools/)
63 - [sklearn](https://scikit-learn.org/stable/)
64 - [MASH](https://github.com/marbl/Mash/releases) (2.3)
65 - [BEDTools](https://bedtools.readthedocs.io/en/latest/) (2.26.0)
66 - [MUmmer](https://github.com/mummer4/mummer) (4.0.0)
67 - [BBTools](https://jgi.doe.gov/data-and-tools/software-tools/bbtools/bb-tools-user-guide/bbmap-guide/) (38.94)
68 - [SKESA](https://github.com/ncbi/SKESA) (2.5.0) [Only required if starting from raw reads]
69
70 ---
71 ## Installing CSP2
72 CSP2 can be installed by cloning the GitHub repo and configuring the [nextflow.config](nextflow.config) and [profiles.config](conf/profiles.config) to suit your needs
73
74 ```
75 git clone https://github.com/CFSAN-Biostatistics/CSP2.git
76 ```
77
78 ## Tips for configuring CSP2
79 CSP2 options can be specified on the command line, or through the Nextflow configuration files detailed in the next section. Feel free to skip this section if you're familiar with editing Nextflow configuration files.
80
81 There are two main configuration files associated with CSP2:
82
83 - The profiles.config file is where you add custom information about your computing environment, but you can also set parameters here as well. An example configuration setup (slurmHPC) is provided as a model.
84
85 - In this example profile, access to the required programs relies on the loading of modules. **However**, there is no need to specify a module for Python, MUMmer, SKESA, bedtools, or MASH if those programs are already in your path.
86
87 ```
88 profiles {
89 standard {
90 process.executor = 'local'
91 params.cores = 1
92 params.python_module = ""
93 params.mummer_module = ""
94 params.skesa_module = ""
95 params.bedtools_module = ""
96 params.mash_module = ""
97 params.bbtools_module = ""
98 }
99 slurmHPC {
100 process.executor = 'slurm'
101 params.cores = 20
102 params.python_module = "/nfs/software/modules/python/3.8.1"
103 params.mummer_module = "/nfs/software/modules/mummer/4.0.0"
104 params.skesa_module = "/nfs/software/modules/skesa/2.5.0"
105 params.bedtools_module = "/nfs/sw/Modules/bedtools"
106 params.bbtools_module = "/nfs/software/modules/bbtools/38.94"
107 params.mash_module = "/nfs/software/modules/mash/2.3"
108 params.trim_name = "_contigs_skesa"
109 }
110 }
111 ```
112 - If you plan to run CSP2 locally, be sure to edit *params.cores* in the standard profile to match the available cores on your system
113 - If you add your own profile, be sure to note it on the command line (one hypen)
114 ```
115 nextflow run CSP2.nf -profile myNewProfile <args>
116 ```
117
118 - The nextflow.config file is where you can change other aspects of the CSP2 run, including data location, QC parameters, and all the options listed below:
119
120 **Options with defaults include**:
121
122 | Parameter | Description | Default Value |
123 |------------------|------------------------------------------------------------------------------------------------------------|-------------------------------------------|
124 | --outroot | Base directory to create output folder | $CWD |
125 | --out | Name of the output folder to create (must not exist) | CSP2_(java.util.Date().getTime()) |
126 | --forward | Full file extension for forward/left reads of query | _1.fastq.gz |
127 | --reverse | Full file extension for reverse/right reads of reference | _2.fastq.gz |
128 | --ref_forward | Full file extension for forward/left reads of reference | _1.fastq.gz |
129 | --ref_reverse | Full file extension for reverse/right reads of reference | _2.fastq.gz |
130 | --readext | Extension for single-end reads for query | fastq.gz |
131 | --ref_readext | Extension for single-end reads for reference | fastq.gz |
132 | --min_cov | Do not analyze queries that cover less than <min_cov>% of the reference assembly | 85 |
133 | --min_iden | Only consider alignments where the percent identity is at least <min_iden>% | 99 |
134 | --min_len | Only consider alignments that span at least <min_len>bp | 500 |
135 | --dwin | A comma-separated list of windows to check SNP densities | 1000,125,15 |
136 | --wsnps | The maximum number of SNPs allowed in the corresponding window from --dwin | 3,2,1 |
137 | --query_edge | Only consider SNPs that occur within <query_edge>bp of the end of a query contig | 150 |
138 | --ref_edge | Only consider SNPs that occur within <query_edge>bp of the end of a reference contig | 150 |
139 | --n_ref | The number of reference isolates to consider (only applied if CSP2 is choosing references) | 1 |
140 | --rescue | If running in SNP Pipeline mode, rescue edge-filtered SNPs that are not edge filtered in 1+ query | Unset (Do not rescue) |
141
142 **Options without defaults include**:
143 | Parameter | Description |
144 |------------------------|------------------------------------------------------------------------------------------------------------------------------------|
145 | --reads | Location of query read data (Path to directory, or path to file with multiple directories) |
146 | --fasta | Location of query assembly data (Path to directory containing FASTAs, path to FASTA, file with list of multiple FASTA paths) |
147 | --ref_reads | Location of reference read data (Path to directory, or path to file with multiple directories) |
148 | --ref_fasta | Location of reference assembly data (Path to directory containing FASTAs, path to FASTA, path to multiple FASTAs) |
149 | --python_module | Name of Python module if 'module load python' statement is required. |
150 | --mummer_module | Name of MUmmer module if 'module load mummer' statement is required. |
151 | --skesa_module | Name of SKESA module if 'module load skesa' statement is required. |
152 | --bedtools_module | Name of BEDTools module if 'module load bedtools' statement is required. |
153 | --mash_module | Name of MASH module if 'module load mash' statement is required. |
154 | --trim_name | A string in assembly file names that you want to remove from sample IDs (e.g., _contigs_skesa) |
155
156 ---
157
158 ## Examples
159
160 The [CSP2 Test Data repo](https://github.com/CFSAN-Biostatistics/CSP2_TestData) contains small test datasets to ensure things are running as expected. Here are a few examples of how you can use CSP2 in screening mode or in SNP pipeline mode.
161
162 ## Screening Mode (Example)
163 *Situation*: As part of a long-term microbiology experiment, you perform weekly WGS sequencing on isolates as they evolve under different selective conditions. As results, your DNA sequencing facility returns raw WGS reads and assembled genomes.
164
165 During Week 42, analyses start detecting high numbers of mutations, and assembly-based results are not concordant with read-based results. You suspect that either the reads or assembly you were given may be from their lab control strain, but you want to check first.
166
167 **The data**:
168 - Read data
169 - Week_42_Reads_1.fq.gz; Week_42_Reads_2.fq.gz
170 - Assembled data:
171 - Week_42_Assembly.fa
172 - Lab_Control.fasta
173
174 In this case, we want to use *--runmode screen*, because we want to explicitly check if either dataset matches the reference strain.
175
176 - **Note**: By default, CSP2 expects read data as zipped fastqs (fastq.gz), with paired-end reads denoted as _1.fastq.gz and _2.fastg.gz. These settings can be changed:
177 - Permanently in the nextflow.config file
178 - Situationally in profiles.config
179 - Directly on the command line, as in the example below:
180 - **Query Reads**: *--readext*; *--forward*; *--reverse*
181 - **Reference Reads**: *--ref_readext*; *--ref_forward*; *--ref_reverse*
182
183 To run this example locally, where *Nextflow*, *SKESA*, *MUMmer*, *Python*, and *BEDTools* are installed on your path, run:
184
185 ```
186 nextflow run CSP2.nf --out Test_Output/Contamination_Screen --runmode screen --ref_fasta assets/Screen/Assembly/Lab_Control.fasta --fasta assets/Screen/Assembly/Week_42_Assembly.fasta --reads assets/Screen/Reads --forward _1.fq.gz --reverse _2.fq.gz --readext fq.gz
187 ```
188
189 ```
190 nextflow run CSP2.nf // Run CSP2
191 --out Test_Output/Contamination_Screen // Save results to ./Test_Output/Contamination_Screen
192 --runmode screen // Compare each query to the reference
193 --ref_fasta assets/Screen/Assembly/Lab_Control.fasta // Compare all queries to this reference
194 --fasta assets/Screen/Assembly/Week_42_Assembly.fasta // Include this assembly as a query
195 --reads assets/Screen/Reads // Include any read datasets from this directory as queries
196 --forward _1.fq.gz // Forward reads don't match the default '_1.fastq.gz'
197 --reverse _2.fq.gz // Reverse reads don't match the default '_2.fastq.gz'
198 --readext fq.gz // Reads don't match the default 'fastq.gz'
199 ```
200
201 If you're running on an HPC and you need to load modules, you could include your custom profile:
202
203 ```
204 # Load Nextflow module if necessary
205 module load nextflow
206
207 nextflow run CSP2.nf -profile slurmHPC --out Test_Output/Contamination_Screen --runmode screen --ref_fasta assets/Screen/Assembly/Lab_Control.fasta --fasta assets/Screen/Assembly/Week_42_Assembly.fasta --reads assets/Screen/Reads --forward _1.fq.gz --reverse _2.fq.gz --readext fq.gz
208 ```
209 ```
210 nextflow run CSP2.nf // Run CSP2
211 -profile slurmHPC // Choose run profile (**note single hyphen**)
212 --out Test_Output/Contamination_Screen // Save results to ./Test_Output/Contamination_Screen
213 --runmode screen // Compare each query to the reference
214 --ref_fasta assets/Screen/Assembly/Lab_Control.fasta // Compare all queries to this reference
215 --fasta assets/Screen/Assembly/Week_42_Assembly.fasta // Include this assembly as a query
216 --reads assets/Screen/Reads // Include any read datasets from this directory as queries
217 --forward _1.fq.gz // Forward reads don't match the default '_1.fastq.gz'
218 --reverse _2.fq.gz // Reverse reads don't match the default '_2.fastq.gz'
219 --readext fq.gz // Reads don't match the default 'fastq.gz'
220 ```
221
222 ### Output
223
224 If all went well, you should see something like this:
225 <p align="center">
226 <img src="img/Screen_Run.jpg" alt="Nextflow print output for screening run"/>
227 </p>
228
229 From top to bottom, we can see that CSP2:
230 - Found the paired-end reads
231 - Assembled them using *SKESA*
232 - Saved an assembly log
233 - Aligned both queries against the reference using MUMmer
234 - Ran the screening script to generate the output table
235
236 Let's take a look to see what was generated:
237 ```
238 ls Test_Output/Contamination_Screen
239
240 Assemblies/
241 Isolate_Data.tsv
242 logs/
243 MUMmer_Output/
244 Raw_MUMmer_Summary.tsv
245 Screening_Results.tsv
246 snpdiffs/
247 ```
248 -**Note**: This output is available to inspect in [assets/Screen/Output](assets/Screen/Output)
249
250 **Directories**
251
252 - *Assemblies*: Directory where any *SKESA* assemblies are stored
253 - *logs*: Directory where run logs are stored
254 - *MUMmer_Output*: Directory where raw MUMmer .snps, .report, .1coords are stored
255 - *snpdiffs*: Directory where .snpdiffs files are stored
256
257 **Files**
258 - *Isolate_Data.tsv*: A TSV file with basic FASTA stats for each isolate
259
260 | Isolate_ID | Isolate_Type | Assembly_Path | Contig_Count | Assembly_Bases | N50 | N90 | L50 | L90 | SHA256 |
261 |------------------|--------------|-----------------------------------------------------------------------------------------------------------------------------------------|--------------|----------------|-------|-------|-----|-----|------------------------------------------------------------------|
262 | Lab_Control | Reference | /flash/storage/scratch/Robert.Literman/NextFlow/CSP2/assets/Screen/Assembly/Lab_Control.fasta | 254 | 4584986 | 41813 | 13112 | 38 | 110 | aae3a07d055bff2fa66127ca77cae35dd5cce5cc42dafea481787d4010c7dbef |
263 | Week_42_Reads | Query | /flash/storage/scratch/Robert.Literman/NextFlow/CSP2/Test_Output/241024_Test_Output/Contamination_Screen/Assemblies/Week_42_Reads.fasta | 372 | 4561747 | 24332 | 7894 | 56 | 171 | 28eca0ecf14fcc1166ae7236c849acc08ad040cd011fc4331ba124db73601009 |
264 | Week_42_Assembly | Query | /flash/storage/scratch/Robert.Literman/NextFlow/CSP2/assets/Screen/Assembly/Week_42_Assembly.fasta | 747 | 4473771 | 12782 | 3438 | 105 | 360 | 85c216de6e1bb9ccf76e7e5d64931884375d66e765f4c4fe726f3be15eb91563 |
265
266
267 - *Screening_Results_PassQC.tsv*: A TSV file with screening results for all queries that made it through QC
268
269 | Query_ID | Reference_ID | Screen_Category | CSP2_Screen_SNPs | Query_Percent_Aligned | Reference_Percent_Aligned | Query_Contigs | Query_Bases | Reference_Contigs | Reference_Bases | Raw_SNPs | Purged_Length | Purged_Identity | Purged_LengthIdentity | Purged_Invalid | Purged_Indel | Purged_Duplicate | Purged_Het | Purged_Density | Filtered_Query_Edge | Filtered_Ref_Edge | Filtered_Both_Edge | Kmer_Similarity | Shared_Kmers | Query_Unique_Kmers | Reference_Unique_Kmers | MUMmer_gSNPs | MUMmer_gIndels |
270 |------------------|--------------|-----------------|------------------|-----------------------|---------------------------|---------------|-------------|-------------------|-----------------|----------|---------------|-----------------|-----------------------|----------------|--------------|------------------|------------|----------------|---------------------|-------------------|--------------------|-----------------|--------------|--------------------|------------------------|--------------|----------------|
271 | Week_42_Assembly | Lab_Control | Pass | 1 | 99.1 | 96.7 | 747 | 4473771 | 254 | 4584986 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 97.32 | 4449973 | 746 | 121740 | 1 | 0 |
272 | Week_42_Reads | Lab_Control | Pass | 50 | 99.15 | 98.65 | 372 | 4561747 | 254 | 4584986 | 61 | 0 | 0 | 12 | 0 | 2 | 0 | 0 | 0 | 1 | 1 | 0 | 98.51 | 4526310 | 23072 | 45403 | 58 | 1 |
273
274 - **Columns**
275
276 | Column | Contents |
277 |---------------------------|----------------------------------------------------------|
278 | Query_ID | Name of query isolate |
279 | Reference_ID | Name of reference isolate |
280 | Screen_Category | Pass/Fail |
281 | CSP2_Screen_SNPs | SNP distance between query and reference |
282 | Query_Percent_Aligned | Percent of query genome aligned to reference |
283 | Reference_Percent_Aligned | Percent of reference genome aligned to query |
284 | Query_Contigs | Contigs in query assembly |
285 | Query_Bases | Bases in query assembly |
286 | Reference_Contigs | Contigs in reference assembly |
287 | Reference_Bases | Bases in reference assembly |
288 | Raw_SNPs | Raw SNPs prior to filtering |
289 | Purged_Length | SNPs purged due to contigs less than min_len |
290 | Purged_Identity | SNPs purged due to percent identity less than min_iden |
291 | Purged_LengthIdentity | SNPs purged to short + low quality alignments |
292 | Purged_Invalid | SNPs purged due to non-ACTG |
293 | Purged_Indel | SNPs purged as indels |
294 | Purged_Duplicate | SNPs purged as duplicate alignments |
295 | Purged_Het | SNPs purged with heterozygous signal |
296 | Purged_Density | SNPs purged for having more than wsnps in dwin |
297 | Filtered_Query_Edge | SNPs purged for being too close to query contig edge |
298 | Filtered_Ref_Edge | SNPs purged for being too close to reference contig edge |
299 | Filtered_Both_Edge | SNPs purged for being too close to both contig edges |
300 | Kmer_Similarity | Kmer similarity between query and reference |
301 | Shared_Kmers | Number of shared kmers between query and reference |
302 | Query_Unique_Kmers | Unique kmer count for query |
303 | Reference_Unique_Kmers | Unique kmer count for reference |
304 | MUMmer_gSNPs | Number of MUMmer gSNPs |
305 | MUMmer_gIndels | Number of MUMmer gIndels |
306
307 - SNPDiffs Files
308 - The .snpdiffs files generated by CSP2 have three main components:
309
310 1. **Header**: The header row of a .snpdiffs file contains all the same information as the screening results TSV.
311 - To peek at a .snpdiffs header, try:
312 ```
313 head -1 (QUERY)__vs__(REF).snpdiffs | tr "\t" "\n"
314 ```
315
316 2. **BED File**: Assuming sufficient overlap in assemblies, the next section will be a BED file of all the 1-to-1 overlaps between the query and reference. This section is denoted by '##\t'
317 ```
318 grep "##" Week_42_Reads__vs__Lab_Control.snpdiffs | head -10
319 ```
320
321
322 | | Reference_Contig | Reference_Align_Start | Reference_Align_End | Reference_Contig_Length | Reference_Contig_Aligned | Query_Contig | Query_Align_Start | Query_Align_End | Query_Contig_Length | Query_Contig_Aligned | Percent_Identity |
323 |----|-------------------------|-----------------------|---------------------|-------------------------|--------------------------|--------------------|-------------------|-----------------|---------------------|----------------------|------------------|
324 | ## | SRR16119110_100_32.8137 | 0 | 6519 | 36207 | 6519 | Contig_277_13.5424 | 4804 | 11323 | 11323 | 6519 | 100 |
325 | ## | SRR16119110_100_32.8137 | 6586 | 7524 | 36207 | 938 | Contig_73_15.3932 | 0 | 938 | 938 | 938 | 100 |
326 | ## | SRR16119110_100_32.8137 | 7828 | 36104 | 36207 | 28276 | Contig_8_13.8951 | 0 | 28276 | 28276 | 28276 | 100 |
327 | ## | SRR16119110_102_54.229 | 225 | 12863 | 12876 | 12638 | Contig_232_18.7715 | 0 | 12638 | 12638 | 12638 | 100 |
328 | ## | SRR16119110_103_30.5388 | 143 | 6585 | 6585 | 6442 | Contig_355_13.8058 | 6300 | 12742 | 12742 | 6442 | 100 |
329 | ## | SRR16119110_104_40.3209 | 0 | 12590 | 19820 | 12590 | Contig_153_15.5804 | 8748 | 21338 | 21338 | 12590 | 100 |
330 | ## | SRR16119110_104_40.3209 | 12664 | 19615 | 19820 | 6951 | Contig_35_15.127 | 0 | 6951 | 6951 | 6951 | 100 |
331 | ## | SRR16119110_105_32.917 | 33 | 11125 | 11125 | 11092 | Contig_57_14.4746 | 37635 | 48727 | 48727 | 11092 | 100 |
332 | ## | SRR16119110_106_22.2235 | 0 | 2959 | 2959 | 2959 | Contig_166_11.9238 | 0 | 2959 | 3512 | 2959 | 100 |
333 | ## | SRR16119110_107_44.0162 | 139 | 12900 | 204844 | 12761 | Contig_174_15.2027 | 0 | 12761 | 12761 | 12761 | 100 |
334
335 3. **SNP Data**: Finally, the third section of a .snpdiffs file contains all the data on SNPs that passed and failed QC. Only SNPs with the Category "SNP" are counted in the final distance.
336 ```
337 grep -v "#" Week_42_Reads__vs__Lab_Control.snpdiffs | head -10
338 ```
339
340 | Reference_Contig | Reference_Loc_Start | Reference_Loc_End | Query_Contig | Query_Loc_Start | Query_Loc_End | Ref_Loc_ID | Query_Loc_ID | Reference_Alignment_Start | Reference_Alignment_End | Query_Alignment_Start | Query_Alignment_End | Reference_Base | Query_Base | Distance_to_Reference_Contig_Edge | Distance_to_Query_Contig_Edge | Reference_Aligned | Query_Aligned | Direction | Percent_Identity | Category |
341 |-------------------------|---------------------|-------------------|--------------------|-----------------|---------------|--------------------------------|--------------------------|---------------------------|-------------------------|-----------------------|---------------------|----------------|------------|-----------------------------------|-------------------------------|-------------------|---------------|-----------|------------------|----------|
342 | SRR16119110_107_44.0162 | 79353 | 79354 | Contig_138_15.1806 | 38294 | 38295 | SRR16119110_107_44.0162/79354 | Contig_138_15.1806/38295 | 41059 | 116706 | 0 | 75647 | A | T | 79354 | 37352 | 75647 | 75647 | 1 | 99.99 | SNP |
343 | SRR16119110_107_44.0162 | 96797 | 96798 | Contig_138_15.1806 | 55738 | 55739 | SRR16119110_107_44.0162/96798 | Contig_138_15.1806/55739 | 41059 | 116706 | 0 | 75647 | C | T | 96798 | 19908 | 75647 | 75647 | 1 | 99.99 | SNP |
344 | SRR16119110_107_44.0162 | 105821 | 105822 | Contig_138_15.1806 | 64762 | 64763 | SRR16119110_107_44.0162/105822 | Contig_138_15.1806/64763 | 41059 | 116706 | 0 | 75647 | G | A | 99022 | 10884 | 75647 | 75647 | 1 | 99.99 | SNP |
345 | SRR16119110_107_44.0162 | 170608 | 170609 | Contig_93_15.202 | 49811 | 49812 | SRR16119110_107_44.0162/170609 | Contig_93_15.202/49812 | 120797 | 204844 | 0 | 84045 | G | A | 34235 | 47655 | 84047 | 84045 | 1 | 99.99 | SNP |
346 | SRR16119110_107_44.0162 | 183055 | 183056 | Contig_93_15.202 | 62258 | 62259 | SRR16119110_107_44.0162/183056 | Contig_93_15.202/62259 | 120797 | 204844 | 0 | 84045 | C | T | 21788 | 35208 | 84047 | 84045 | 1 | 99.99 | SNP |
347 | SRR16119110_135_26.5159 | 2370 | 2371 | Contig_215_11.8373 | 2370 | 2371 | SRR16119110_135_26.5159/2371 | Contig_215_11.8373/2371 | 0 | 2511 | 0 | 2511 | G | A | 140 | 789 | 2511 | 2511 | 1 | 99.96 | SNP |
348 | SRR16119110_150_30.823 | 36612 | 36613 | Contig_168_14.004 | 36883 | 36884 | SRR16119110_150_30.823/36613 | Contig_168_14.004/36884 | 0 | 49916 | 271 | 50187 | C | T | 13485 | 13303 | 49916 | 49916 | 1 | 99.99 | SNP |
349 | SRR16119110_153_30.6937 | 3685 | 3686 | Contig_233_12.8229 | 3685 | 3686 | SRR16119110_153_30.6937/3686 | Contig_233_12.8229/3686 | 0 | 7254 | 0 | 7254 | G | A | 3568 | 3594 | 7254 | 7254 | 1 | 99.99 | SNP |
350 | SRR16119110_156_32.1591 | 12089 | 12090 | Contig_281_12.9967 | 13905 | 13906 | SRR16119110_156_32.1591/12090 | Contig_281_12.9967/13906 | 0 | 24120 | 1816 | 25936 | A | T | 12030 | 13591 | 24120 | 24120 | 1 | 99.99 | SNP |
351 | SRR16119110_156_32.1591 | 19366 | 19367 | Contig_281_12.9967 | 21182 | 21183 | SRR16119110_156_32.1591/19367 | Contig_281_12.9967/21183 | 0 | 24120 | 1816 | 25936 | A | G | 4753 | 6314 | 24120 | 24120 | 1 | 99.99 | SNP |
352
353 **Conclusions**
354
355 Based on these results, the assembly provided by the DNA sequencing facility had only 1 SNP relative to the lab control strain, but the read data was over 40 SNPs away, suggesting that the wrong assembly was packaged with the read data.
356
357 ## SNP Pipeline Mode (Example)
358 *Situation*: You dug 10 soil samples and isolated a single microbe from each. You want to check the relatedness of the cultured isolates.
359
360 The data:
361 - Assemblies from soil isolates (Sample_(A-O).fasta)
362
363 In this case, we want to use *--runmode snp*, because we want to calculate the pairwise distances between all isolates and generate alignments. We will let RefChooser choose the best reference genome.
364
365 ```
366 nextflow run CSP2.nf --out Test_Output/Soil_Analysis --runmode snp --fasta assets/SNP/
367 ```
368 ```
369 nextflow run CSP2.nf // Run CSP2
370 --out Test_Output/Soil_Analysis // Save results to ./Test_Output/Soil_Analysis
371 --runmode snp // Compare all queries to each other
372 --fasta assets/SNP // Gather query assemblies from this directory (you can also point to a text file containing multiple FASTA paths)
373 ```
374
375 ### Output
376
377 If all went well, you should see something like this:
378 <p align="center">
379 <img src="img/SNP.jpg" alt="Nextflow print output for SNP Pipeline run"/>
380 </p>
381
382 From top to bottom, we can see that CSP2:
383 - Skipped assembly steps, as queries were already assembled
384 - Ran *MASH* to choose a suitable reference genome
385 - Aligned each query to the reference using MUMmer
386 - Ran the SNP pipeline workflow on the .snpdiffs files to generate alignments and SNP distance data
387
388 Let's take a look to see what was generated:
389 ```
390 ls Test_Output/Soil_Analysis
391
392 Isolate_Data.tsv
393 logs/
394 MUMmer_Output/
395 Raw_MUMmer_Summary.tsv
396 SNP_Analysis/
397 snpdiffs/
398 ```
399 -**Note**: This output is available to inspect in [assets/SNP/Output](assets/SNP/Output)
400
401 **Directories**
402 - *logs*: Directory where run logs are stored
403 - *MUMmer_Output*: Directory where raw MUMmer .snps, .report, .1coords are stored
404 - *snpdiffs*: Directory where .snpdiffs files are stored
405 - *SNP_Analysis*: Directory where SNP pipeline results are stored
406
407 **Files**
408
409 - The log file for any SNP pipeline run contains lots of useful information, including runtimes and output locations.
410
411 ```
412 cat logs/SNP_Sample_A.log
413
414 CSP2 SNP Pipeline Analysis
415 Reference Isolate: Sample_A
416 2024-10-17 11:37:44
417 -------------------------------------------------------
418
419 Reading in SNPDiffs files...Done!
420 - Read in 14 SNPdiffs files
421 -------------------------------------------------------
422
423 - Not performing SNP edge rescuing (any SNPs found within 150bp of a query contig edge will be purged)...
424 -------------------------------------------------------
425
426 Screening all queries against reference...Done!
427 - Reference screening data saved to /flash/storage/scratch/Robert.Literman/NextFlow/Test_CSP2/241024_Test_Output/Soil_Analysis/SNP_Analysis/Sample_A/Reference_Screening.tsv
428 - Of 14 comparisons, 14 covered at least 85.0% of the reference genome after removing poor alignments
429 -------------------------------------------------------
430
431 Compiling SNPs across 14 samples...
432 - 220 total SNPs detected across all samples...
433 - 167 unique SNPs passed QC filtering in at least one sample...
434 - 1 unique SNPs were within 150bp of a reference contig end and were not considered in any query...
435 - Skipping edge resucing...
436 - 53 unique SNPs were purged in all queries they were found in, and were not considered in the final dataset...
437 - Processed coverage information...
438 - SNP coverage information: /flash/storage/scratch/Robert.Literman/NextFlow/Test_CSP2/241024_Test_Output/Soil_Analysis/SNP_Analysis/Sample_A/Locus_Categories.tsv
439 - Query coverage information: /flash/storage/scratch/Robert.Literman/NextFlow/Test_CSP2/241024_Test_Output/Soil_Analysis/SNP_Analysis/Sample_A/Query_Coverage.tsv
440 -------------------------------------------------------
441
442 Processing alignment data...Done!
443 - Saved alignment of 167 SNPs to /flash/storage/scratch/Robert.Literman/NextFlow/Test_CSP2/241024_Test_Output/Soil_Analysis/SNP_Analysis/Sample_A/snpma.fasta
444 - Saved ordered loc list to /flash/storage/scratch/Robert.Literman/NextFlow/Test_CSP2/241024_Test_Output/Soil_Analysis/SNP_Analysis/Sample_A/snplist.txt
445 - Preserving SNPs with at most 50.0% missing data...
446 - Of 167 SNPs, 167 SNPs pass the 50.0% missing data threshold...
447 - Saved preserved alignment to /flash/storage/scratch/Robert.Literman/NextFlow/Test_CSP2/241024_Test_Output/Soil_Analysis/SNP_Analysis/Sample_A/snpma_preserved.fasta
448 - Saved preserved ordered loc list to /flash/storage/scratch/Robert.Literman/NextFlow/Test_CSP2/241024_Test_Output/Soil_Analysis/SNP_Analysis/Sample_A/snplist_preserved.txt
449 -------------------------------------------------------
450
451 Processing pairwise comparisons files...Done!
452 - Saved raw pairwise distances to /flash/storage/scratch/Robert.Literman/NextFlow/Test_CSP2/241024_Test_Output/Soil_Analysis/SNP_Analysis/Sample_A/snp_distance_pairwise.tsv
453 - Saved raw pairwise matrix to /flash/storage/scratch/Robert.Literman/NextFlow/Test_CSP2/241024_Test_Output/Soil_Analysis/SNP_Analysis/Sample_A/snp_distance_matrix.tsv
454 - Saved preserved pairwise distances to /flash/storage/scratch/Robert.Literman/NextFlow/Test_CSP2/241024_Test_Output/Soil_Analysis/SNP_Analysis/Sample_A/snp_distance_pairwise_preserved.tsv
455 - Saved preserved pairwise matrix to /flash/storage/scratch/Robert.Literman/NextFlow/Test_CSP2/241024_Test_Output/Soil_Analysis/SNP_Analysis/Sample_A/snp_distance_matrix_preserved.tsv
456 Total Time: 1.45 seconds
457 -------------------------------------------------------
458 ```
459
460 - For each reference sample there will be a folder called SNP_Analysis/REFERENCE_ID, containing:
461 - *CSP2_SNP_Pipeline.log*: Log file for SNP Pipeline run
462 - *Locus_Categories.tsv*: A TSV file containing information about each SNP and how many isolates had missing or purged data
463 - *Query_Coverage.tsv*: A TSV file containing information about each query and how many sites had missing or purged data
464 - *Reference_Screening.tsv*: Output from CSP2 Screening mode
465 - *snpdiffs.txt*: File containing a list of all *snpdiffs* files used in the analysis
466 - *snp_distance_matrix.tsv*: SNP distances between isolates in matrix format
467 - *snp_distance_pairwise.tsv*: SNP distances between isolates in column format
468 - *snplist.txt*: List of reference loc IDs (Reference_Contig/End_Position)
469 - *snpma.fasta*: Alignment of SNP data
470
471 - If using the --max_missing argument, CSP2 will also output files where SNPs with too many missing data are removed:
472
473 - *snp_distance_matrix_preserved.tsv*
474 - *snp_distance_pairwise_preserved.tsv*
475 - *snplist_preserved.txt*
476 - *snpma_preserved.fasta*
477
478 To see if the reference genome has an impact on SNP distance estimation, you can test one or more via *--ref_fasta* / *--ref_reads* or have CSP2 choose *--n_ref* isolates
479 ```
480 nextflow run CSP2.nf --out Test_Output/Soil_Analysis --runmode snp --fasta assets/SNP/ --n_ref 3
481 ```
482 ```
483 nextflow run CSP2.nf // Run CSP2
484 --out Test_Output/Soil_Analysis // Save results to ./Test_Output/Soil_Analysis
485 --runmode snp // Compare all queries to each other
486 --fasta assets/SNP // Gather query assemblies from this directory
487 --n_ref 3 // Choose the top 3 references and run 3 sepearate analyses
488 ```