SciCodePile/SciCode-Domain-Code
DATA1: Domain-Specific Code Dataset Dataset Overview DATA1 is a large-scale domain-specific code dataset focusing on code samples from interdisciplinary fields such as biology, chemistry, materials science, and related areas. The dataset is collected and organized from GitHub repositories, covering 178 different domain topics with over 1.1 billion lines of code. Dataset Statistics Total Datasets: 178 CSV files Total Data Size: ~115 GB Total Lines… See the full description on the dataset page: https://huggingface.co/datasets/SciCodePile/SciCode-Domain-Code.
41.3k
1"keyword","repo_name","file_path","file_extension","file_size","line_count","content","language"
2"Microbiology","microsud/Tools-Microbiome-Analysis","CODE_OF_CONDUCT.md",".md","3359","77","# Contributor Covenant Code of Conduct3 4## Our Pledge5 6In the interest of fostering an open and welcoming environment, we as7contributors and maintainers pledge to making participation in our project and8our community a harassment-free experience for everyone, regardless of age, body9size, disability, ethnicity, sex characteristics, gender identity and expression,10level of experience, education, socio-economic status, nationality, personal11appearance, race, religion, or sexual identity and orientation.12 13## Our Standards14 15Examples of behaviour that contributes to creating a positive environment16include:17 18* Using welcoming and inclusive language19* Being respectful of differing viewpoints and experiences20* Gracefully accepting constructive criticism21* Focusing on what is best for the community22* Showing empathy towards other community members23 24Examples of unacceptable behavior by participants include:25 26* The use of sexualized language or imagery and unwelcome sexual attention or27 advances28* Trolling, insulting/derogatory comments, and personal or political attacks29* Public or private harassment30* Publishing others' private information, such as a physical or electronic31 address, without explicit permission32* Other conduct which could reasonably be considered inappropriate in a33 professional setting34 35## Our Responsibilities36 37Project maintainers are responsible for clarifying the standards of acceptable38behavior and are expected to take appropriate and fair corrective action in39response to any instances of unacceptable behavior.40 41Project maintainers have the right and responsibility to remove, edit, or42reject comments, commits, code, wiki edits, issues, and other contributions43that are not aligned to this Code of Conduct, or to ban temporarily or44permanently any contributor for other behaviors that they deem inappropriate,45threatening, offensive, or harmful.46 47## Scope48 49This Code of Conduct applies both within project spaces and in public spaces50when an individual is representing the project or its community. Examples of51representing a project or community include using an official project e-mail52address, posting via an official social media account, or acting as an appointed53representative at an online or offline event. Representation of a project may be54further defined and clarified by project maintainers.55 56## Enforcement57 58Instances of abusive, harassing, or otherwise unacceptable behavior may be59reported by contacting the project team at sudarshanshetty9@gmail.com. All60complaints will be reviewed and investigated and will result in a response that61is deemed necessary and appropriate to the circumstances. The project team is62obligated to maintain confidentiality with regard to the reporter of an incident.63Further details of specific enforcement policies may be posted separately.64 65Project maintainers who do not follow or enforce the Code of Conduct in good66faith may face temporary or permanent repercussions as determined by other67members of the project's leadership.68 69## Attribution70 71This Code of Conduct is adapted from the [Contributor Covenant][homepage], version 1.4,72available at https://www.contributor-covenant.org/version/1/4/code-of-conduct.html73 74[homepage]: https://www.contributor-covenant.org75 76For answers to common questions about this code of conduct, see77https://www.contributor-covenant.org/faq78","Markdown"
79"Microbiology","microsud/Tools-Microbiome-Analysis","setup_microbiome_analysis.R",".R","2908","87","80############################################################################################################################################81# #82# START OF CODE #83# #84############################################################################################################################################85 86# These are some R pacakges I commonly use for analysis. Everytime I update R to latest version, I run this installation to in setup my environment. 87 88# Copy from here89 90setup_microbiome_analysis <- function(){91 92 .packages = c(""ape"", 93 ""gridExtra"", 94 ""picante"", 95 ""data.table"", 96 ""RColorBrewer"", 97 ""DT"", 98 ""reshape"", 99 ""reshape2"", 100 ""magrittr"", 101 ""markdown"",102 ""ggpubr"", 103 ""tibble"", 104 ""pheatmap"", 105 ""dplyr"", 106 ""viridis"", 107 ""devtools"", 108 ""rmdformats"",109 ""intergraph"",110 ""network"",111 ""igraph"",112 ""ggplot2"", 113 ""gridExtra"", 114 ""knitr"", 115 ""vegan"", 116 ""plyr"", 117 ""dplyr"",118 ""ggrepel"", 119 ""ggnetwork"", 120 ""ade4"", 121 ""rmarkdown"",122 ""formatR"",123 ""caTools"",124 ""GGally"")125 126 .bioc_packages <- c(""phyloseq"",127 ""microbiome"", 128 ""phangorn"", 129 ""genefilter"")130 if (!requireNamespace(""BiocManager"", quietly = TRUE))131 install.packages(""BiocManager"")132 133 # Install CRAN packages (if not already installed)134 .inst <- .packages %in% installed.packages()135 if(length(.packages[!.inst]) > 0) install.packages(.packages[!.inst])136 137 .inst <- .bioc_packages %in% installed.packages()138 if(any(!.inst)) {139 BiocManager::install(.bioc_packages[!.inst])140 }141 142 143 if (!(""microbiomeutilities"" %in% installed.packages())) {144 devtools::install_github(""microsud/microbiomeutilities"")145 }146 147 if (!(""SpiecEasi"" %in% installed.packages())) {148 devtools::install_github(""zdk123/SpiecEasi"")149 }150 151 if (!(""ggnet"" %in% installed.packages())) {152 devtools::install_github(""briatte/ggnet"")153 }154 155 message(""If there was no error then you are ready to do microbiome data analysis"")156 157}158 159setup_microbiome_analysis()160 161# Copy until here previous line!162 163####################################################END OF CODE###########################################################################164 165","R"
166"Microbiology","lskatz/lyve-SET","plugins/correct_vcf_filter.py",".py","2241","71","#!/usr/local/bin/python167# Khushbu Patel | 05/15/2018168# Corrects the filter and Calculates number of bases falsely assigned incorrect filter when the bases pass the consensus and coverage. Prints out a new vcf with filters corrected. 169# Python 3.4170# Usage: ./correct_vcf_filter.py inputfile.vcf171# Output will be another vcf file, named correct_vcf.vcf; This script will not overwrite the original VCF172 173import sys174import os175import subprocess176 177infile= sys.argv[1] # Takes the file name as a command line argument178f=open(infile,""r"")179 180basename = os.path.basename(infile)181outfile_vcf = basename + ""_corrected.vcf""182outfile_txt = basename + ""_corrected.txt""183temp = []184count = 0185out = []186ID = ''187old_filter = 0188 189 190with open(outfile_vcf, 'a') as f1:191 with open(outfile_txt, 'a') as f2:192 f2.write( ""CHROM\tPOS\tNEW_FILTER\tDEPTH\tFREQ\tOLD_FILTER\n"")193 for line in f:194 line =line.rstrip()195 if line.startswith('#'):196 ID = line # Printing the headers197 f1.write(ID)198 f1.write(""\n"")199 200 else:201 array = line.split()202 temp = array[9].split(':')203 204 if(array[6] != ""PASS""):205 if(len(temp)> 6):206 old_filter = array[6] # Stores old filter207 temp[6] = temp[6].replace('%','')208 if(int(temp[3]) >= 20 and float(temp[6]) < 5.0): # if coverage and consensus both meet, change the filter to PASS209 #print(line, ""--Incorrect Filter!"") # Sanity check!210 array[6] = ""PASS""211 count += 1 # Count the number of sites that have been assigned wrong filter212 out = '\t'.join(array)213 f1.write(out)214 str = array[0]+'\t'+array[1]+'\t'+array[6]+'\t'+temp[3]+'\t'+temp[6]+'\t'+old_filter+'\n'215 f2.write(str)216 217 218 else: # if allele frequency is greater than 5% or coverage does not meet; keep the filter219 out = '\t'.join(array)220 f1.write(out)221 f1.write(""\n"")222 223 else: # else - print everything as is224 out = '\t'.join(array)225 f1.write(out)226 f1.write(""\n"")227 228 else: # else - print everything as is229 out = '\t'.join(array)230 f1.write(out)231 f1.write(""\n"")232 233 234 235print('Number of sites that had been assigned wrong filters %d'%count)236","Python"
237"Microbiology","lskatz/lyve-SET","docs/TESTDATA.md",".md","999","16","To make a test dataset, create a directory with all lowercase with the the files SRA and NUCLEOTIDE. These files are described below; you do not have to create any other files in the test data directory.238 239All directories must start off by containing the files SRA and NUCLEOTIDE which have identifiers for raw reads and for genome assemblies (draft or finished). The first identifier in NUCLEOTIDE will be the reference genome.240 241The final files/directories in each must be the following242asm/ NUCLEOTIDE reads/ reference/ SRA243* asm/ has all the genome assemblies244* NUCLEOTIDE has all assembly identifiers. The reference genome is the first listing245* reads/ has all raw read fastq files246* reference/ contains a single reference genome assembly247* SRA has a listing of all raw read identifiers that should be downloaded248 249CREDITS250* Test lambda data were obtained from the CFSAN SNP pipeline at https://github.com/CFSAN-Biostatistics/snp-pipeline251* Other datasets were created between EDLB and CFSAN252","Markdown"
253"Microbiology","lskatz/lyve-SET","docs/OUTPUT.md",".md","7316","133","Below is a visualization of the workflow with output files.254Then in the next section, a table with the description of all output files.255 256Visualization of output files257=============================258 259```mermaid260flowchart TD261 subgraph LOGDIR [""log directory""]262 LOG[""main log file""]263 subgraph LOGSUBDIR [""SGELK subfolder""]264 LOGS[""launch_set.pl.$$.log: individual log files""]265 JOBS[""qsub.$$.pl: individual jobs to execute""]266 FINISHED[""launch_set.pl.$$.finished: individual jobs that were finished""]267 SUBMITTED[""launch_set.pl.$$.submitted: individual jobs that were submitted""]268 RUNNING[""launch_set.pl.$$.running: individual jobs that are running""]269 end270 end271 272 SET_MANAGE_CREATE[--create] --> |Create all directories such as reads, reference| READSDIR273 SET_MANAGE_READS[--add-reads] --> |Create all directories such as reads, reference| READSDIR274 SET_MANAGE_ASM[--add-assembly] --> |Add an assembly| ASMDIR275 SET_MANAGE_REF[--change-reference] --> |Add a reference genome| REFDIR276 SET_MANAGE --- SET_MANAGE_CREATE277 SET_MANAGE --- SET_MANAGE_READS278 SET_MANAGE --- SET_MANAGE_ASM279 SET_MANAGE --- SET_MANAGE_REF280 subgraph READSDIR [""reads directory""]281 direction LR;282 SEQUENCER[""R1 R2 fastq files""]283 RAW[""Raw reads""]284 CLEANED[""Cleaned reads""]285 286 RAW --> |run_assembly_trimClean.pl| CLEANED287 SEQUENCER --> |shuffleSplitReads.pl| RAW288 289 end290 subgraph REFDIR [""reference directory""]291 direction LR;292 REF[""Reference genome""]293 UNMASKEDBED[""unmaskedRegions.bed""]294 MASKEDBED[""maskedRegions.bed""]295 296 REF --> |findPhages.pl| MASKEDBED297 MASKEDBED --> |invert| UNMASKEDBED298 end299 subgraph ASMDIR [""asm directory""]300 INASM[""input assemblies""]301 end302 LAUNCH_SMALT{{launch_smalt.pl}}303 304 ASMDIR --> |samtools wgsim| READSDIR305 READSDIR --> |launch_smalt.pl| LAUNCH_SMALT306 REFDIR --> |launch_smalt.pl; only accept unmasked regions| LAUNCH_SMALT307 LAUNCH_SMALT --> |make a bam, once per genome| BAMDIR308 subgraph BAMDIR [""bam directory""]309 direction LR;310 BAM[""*.bam""]311 BAMIDX[""*.bam.bai""]312 BAMCLIFFS[""*.bam.cliffs.bed""]313 314 BAM --> |set_findCliffs.pl| BAMCLIFFS315 BAM --> |samtools index| BAMIDX316 end317 LAUNCH_VARSCAN{{launch_varscan.pl}}318 319 BAMDIR --> |launch_varscan.pl \nexclude any regions in *.bam.cliffs.bed\nAccept SNPs at %consenus, X depth, fwd/rev support| LAUNCH_VARSCAN320 REFDIR --> |launch_varscan.pl\ndo not reprocess with unmaskedRegions.bed \nb/c launch_smalt.pl already used it| LAUNCH_VARSCAN321 LAUNCH_VARSCAN --> |make a vcf, once per genome| VCFDIR322 subgraph VCFDIR [""VCF directory""]323 direction LR;324 VCF[""vcf files""]325 VCFIDX[""vcf index files""]326 327 VCF --> |bcftools index| VCFIDX328 end329 VCFDIR --> |set_mergeVcf.sh: combine all Vcfs into \nout.pooled.vcf.gz and out.pooled.snps.vcf.gz| MSADIR330 subgraph MSADIR [""MSA directory""]331 direction LR;332 MERGEVCF{{set_mergeVcf.sh}}333 VCFPOOLED[""out.pooled.vcf.gz""]334 VCFSNPSPOOLED[""out.pooled.snps.vcf.gz""]335 SNPMATRIX[""out.snpmatrix.tsv""]336 FILTEREDMATRIX[""out.filteredMatrix.tsv""]337 FULLALN[""out.aln.fasta""]338 INFORMATIVEALN[""out.informative.fasta""]339 TREE[""out.RAxML_bipartitions; tree.dnd""]340 PAIRWISE[""out.pairwise.tsv""]341 PAIRWISEMATRIX[""out.pairwiseMatrix.tsv""]342 343 MERGEVCF --> VCFPOOLED344 MERGEVCF --> VCFSNPSPOOLED345 VCFPOOLED --> |pooledToMatrix.sh| SNPMATRIX346 SNPMATRIX --> |filterMatrix.pl| FILTEREDMATRIX347 SNPMATRIX --> |matrixToAlignment.pl| FULLALN348 FILTEREDMATRIX --> |matrixToAlignment.pl| INFORMATIVEALN349 FULLALN --> |pairwiseDistances.pl| PAIRWISE350 PAIRWISE --> |pairwiseTo2d.pl| PAIRWISEMATRIX351 INFORMATIVEALN --> |launch_raxml.sh| TREE352 end353```354 355Output files356============357| File | Description | Notes |358|:--------|:---------------|:------|359|project/msa | The multiple sequence alignment directory | Most of the output files you want are here like the multiple sequence alignment and the phylogeny|360|`project/msa/out.pooled.vcf.gz` | The pooled VCF file created from `bcftools merge` | 361|`project/msa/out.pooled.snps.vcf.gz` | SNPs vcf | The same data as `out.pooled.vcf.gz` but filtered to SNPs only. |362|`project/msa/out.pooled.vcf.gz.tbi`, `out.pooled.snps.vcf.gz.tbi` | the tabix index file for each VCF | 363|`project/msa/out.snpmatrix.tsv` | The `bcftools query` output | This file is essentially the main SNP matrix and describes the position and allele for each genome. Each allele is in the genotype (GT) format, as specified in the vcf format specification |364|`project/msa/out.filteredMatrix.tsv` | The filtered `bcftools query` output | After `out.snpmatrix.tsv` is generated, this file describes remaining SNPs after some are filtered out, usually because the `--allowedFlanking` option in `launch_set.pl`, `--allowed` in `filterMatrix.pl`, or similar parameters in other scripts |365|`project/msa/out.aln.fasta` | The output alignment file in fasta format. | Make any changes to this file before running a phylogeny |program. Do not use `out.informative.fasta` to make edits because positions might come and go and therefore you might lose resolution. After any edits, use `removeUninformativeSites.pl` to re-create `out.informative.fasta` |366| `project/msa/out.informative.fasta` | The alignment after removing uninformative columns (ambiguities, invariants, gaps) | Do not make any changes to this file before running a phylogeny. Make the changes in `out.aln.fasta` |367| `project/msa/out.RAxML_bipartitions` | RAxML-generated tree in newick format | 368| `project/msa/tree.dnd` | Symlink to `out.RAxML_bipartitions`| 369| `project/msa/out.pairwise.tsv` | Pairwise distances file | Format: tab-delimited with three columns: genome1, genome2, hqSNP distance |370| `project/msa/out.pairwiseMatrix.tsv` | Pairwise distances matrix | The same data as `out.pairwise.tsv`, but in a 2-d matrix. Generated with `pairwiseTo2d.pl`. |371|project/log| Log files|372|`project/log/launch_set.log` | The main log file |373|project/asm, project/reads || The input assemblies and reads. |374|project/reference | Where the reference fasta file is|375|`project/reference/maskedRegions.bed` | Regions of the reference genome that is masked for analysis. |376|project/reference/maskedRegions | BED-formatted files that describe regions that should be masked in the reference genome.| You may also create your own file that can have any filename with extension `.bed`. This file can describe your manually-chosen regions that should be masked. These regions will be incorporated into `project/reference/maskedRegions.bed`.|377|`project/reference/maskedRegions/phages.bed`| BED-formatted file describing predicted phage sites||378|project/bam| Output bam files are here|379|`project/bam/*.sorted.bam` | Sorted bam files | The query and reference name are encoded in the filename; many times the reference name will just be called ""reference."" |380|`project/bam/*.sorted.bam.bai` | Samtools index file |381|`project/bam/*.sorted.bam.cliffs.bed` | Files describing genome depth cliffs | These are only present if you specified `--mask-cliffs` |382|project/vcf |VCF files|Have the same file format as the `*.sorted.bam` files, so that they can be matched easily when running Lyve-SET. These files are sorted with vcftools and compressed with bgzip.|383|`project/vcf/*.vcf.gz`|VCF files ||384|`project/vcf/*.vcf.gz.tbi`| Tabix index files|385","Markdown"
386"Microbiology","lskatz/lyve-SET","docs/TIPS.md",".md","16327","348","# Tips and Tricks387 388Here are just some tips and tricks that I have used or that others have contributed389 390## Masking a region in your reference genome391 392Yes, there is actually a mechanism to manually mask troublesome regions in the reference genome! Under `project/reference/maskedRegions`, create a file with an extension `.bed`. This file has at least three columns: `contig`, `start`, `stop`. BED is a standard file format and is better described here: https://genome.ucsc.edu/FAQ/FAQformat.html#format1 393 394The Lyve-SET phage-finding tool that uses PHAST actually puts a `phages.bed` file into that directory. In the course of the pipeline, Lyve-SET will use any BED files in that directory to 1) ignore any reads that are mapped entirely in those regions and 2) ignore any SNPs that are found in those regions. In the future, Lyve-SET will also use individualized BED files in the bam directory to mask SNPs found on a per-genome basis.395 396## Using multiple processors on a single-node machine397 398Unfortunately if you are not on a cluster, then Lyve-SET will only work on a single node. Here are some ways of using xargs to speed up certain steps of Lyve-SET. Incorporating these changes or something similar is on my todo list but for now it is easier to post them.399 400### Making pileups401 402In the case where you want to generate all the pileups on one node using SNAP using 24 cores403 404 ls reads/*.fastq.gz | xargs -n 1 -I {} basename {} | xargs -P 24 -n 1 -I {} launch_smalt.pl -f reads/{} -b bam/{}-2010EL-1786.sorted.bam -t tmp -r reference/2010EL-1786.fasta --numcpus 1405 406### Calling SNPs407 408In the case where you have all pileups finished and want to call SNPs on a single node. This example uses 24 cpus. At the end of this example, you will still need to sort the VCF files (`vcf-sort`), compress them with `bgzip`, and index them with `tabix`.409 410 # Call SNPs into vcf file411 ls bam/*.sorted.bam | xargs -n 1 -I {} basename {} .sorted.bam | xargs -P 24 -n 1 -I {} sh -c ""/home/lkatz/bin/Lyve-SET/scripts/launch_varscan.pl bam/{}.sorted.bam --tempdir tmp --reference reference/2010EL-1786.fasta --altfreq 0.75 --coverage 10 > vcf/{}.vcf""412 # sort/compress/index413 cd vcf; 414 ls *.vcf| xargs -I {} -P 24 -n 1 sh -c ""vcf-sort < {} > {}.tmp && mv -v {}.tmp {} && bgzip {} && tabix {}.gz""415 416## SNP interrogation417 418### What are the differences between just two isolates?419 420Sometimes you just want to know the SNPs that define one isolate, or maybe find which SNPs define which clades.421Here is a method to show unique SNPs between two isolates.422You can expand this idea to more isolates or otherwise refine it as needed.423 424```bash425$ cut -f 1,2,5,6 out.snpmatrix.tsv | awk '$3!=$4 && $3!=""N"" && $4!=""N""' | head -n 3426# [1]CHROM [2]POS [5]sample1:GT [6]sample2:GT427gi|9626243|ref|NC_001416.1| 403 A G428gi|9626243|ref|NC_001416.1| 753 A G429```430 431In this example, only the first two sites are shown (`head -n 3`)432and show that at positions 403 and 753, there are differences between sample1 and sample2.433Sites with `N` are ignored with logic like `$3!=""N""`.434Otherwise `awk` is just selecting for rows where the samples differ (`$3!=$4`)435`cut` is used to grab only the position columns 1 and 2 and then to keep only two samples at columns 5 and 6.436 437### What are the differences between two groups of isolates?438 439To do this, we can run `bcftools isec` in combination with `bcftools filter` and `bcftools view`.440In this section, we can use the lambda subset to pretend that there are two clades with clade1 having sample1 and sample2,441and clade2 having samples 3 and 4.442 443Go into the `msa` subfolder with `cd msa`.444 445 # Make files with sample1 and sample2, and then sample3 and sample4446 echo -e ""sample1\nsample2"" > clade1.txt447 echo -e ""sample3\nsample4"" > clade2.txt448 449Use the file that already contains SNPs for all samples and divide it into the two clades450 451 bcftools filter -i 'ALT!=""N""' out.pooled.snps.vcf.gz | \452 bcftools view -S clade1.txt | \453 bgzip -c > clade1.snps.vcf.gz454 bcftools filter -i 'ALT!=""N""' out.pooled.snps.vcf.gz | \455 bcftools view -S clade2.txt | \456 bgzip -c > clade2.snps.vcf.gz457 458Index the new VCF files459 460 tabix clade1.snps.vcf.gz461 tabix clade2.snps.vcf.gz462 463Find the intersection of all SNPs between the two clades464 465 bcftools isec -p isec clade1.snps.vcf.gz clade2.snps.vcf.gz466 467Go into the new `isec` folder and format all the vcf files468 469 cd isec470 for i in *.vcf; do471 bgzip $i472 tabix $i.gz473 done474 475View the notes on which file is which. Typically, `0000.vcf.gz` is going to be sites exclusive to only clade1.476`0001.vcf.gz` is going to be sites exclusive to only clade2.477`0002.vcf.gz` is records in clade1 that are also in clade2.478`0003.vcf.gz` is records in clade2 that are also in clade1.479 480 cat README.txt481 482 This file was produced by vcfisec.483 The command line was: bcftools isec -p isec clade1.snps.vcf.gz clade2.snps.vcf.gz484 485 Using the following file names:486 isec/0000.vcf for records private to clade1.snps.vcf.gz487 isec/0001.vcf for records private to clade2.snps.vcf.gz488 isec/0002.vcf for records from clade1.snps.vcf.gz shared by both clade1.snps.vcf.gz clade2.snps.vcf.gz489 isec/0003.vcf for records from clade2.snps.vcf.gz shared by both clade1.snps.vcf.gz clade2.snps.vcf.gz490 491Merge back into one happy VCF file492 493 bcftools merge 0002.vcf.gz 0003.vcf.gz > merged.vcf494 bgzip merged.vcf495 tabix merged.vcf.gz496 497View it. These are now SNPs that _could_ differentiate clade1 from clade2, but not necessarily 100%.498For example, a SNP might only occur in sample1 but not sample2 even though it's the same clade.499 500 query -H -f '%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n' merged.vcf.gz | column -ts $'\t' | less -S501 502In our example with lambda, with its star phylogeny, however, we do not have SNPs that are 100%.503Hopefully your results differ in this regard :) 504 505 # [1]CHROM [2]POS [3]REF [4]ALT [5]sample1:GT [6]sample2:GT [7]sample3:GT [8]sample4:GT506 gi|9626243|ref|NC_001416.1| 403 G A 1/1 0/0 0/0 0/0507 gi|9626243|ref|NC_001416.1| 550 G A 0/0 0/0 0/0 1/1508 gi|9626243|ref|NC_001416.1| 586 C G 0/0 0/0 0/0 1/1509 gi|9626243|ref|NC_001416.1| 753 G A 1/1 0/0 0/0 0/0510 gi|9626243|ref|NC_001416.1| 1019 C T 0/0 0/0 0/0 1/1511 512### SNP counting513 514How many sites are in the Lyve-SET analysis? Count the number of lines in `out.snpmatrix.tsv`. Subtract 1 for the header.515 516How many SNPs are in the Lyve-SET analysis? Count the number of lines in `out.filteredMatrix.tsv`. Subtract 1 for the517header.518 519I want to count the number of sites as defined as X. Use `out.snpmatrix.tsv` and `filterMatrix.pl` to filter it your way. If you are an advanced user, you can use `bcftools query` on `out.pooled.vcf.gz`.520 521### Why did a SNP fail?522 523To see why any SNP failed, view the vcf files in the VCF directory. Parse them with `bcftools`.524 525#### View the FILTER column526 527 $ bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%FILTER\n' sample1.fastq.gz-reference.vcf.gz | head528 gi|9626243|ref|NC_001416.1| 1 G N DP10;AF0.75529 gi|9626243|ref|NC_001416.1| 2 G N DP10;AF0.75530 gi|9626243|ref|NC_001416.1| 3 G N DP10;AF0.75531 gi|9626243|ref|NC_001416.1| 4 C N DP10;AF0.75532 gi|9626243|ref|NC_001416.1| 5 G N DP10;AF0.75533 gi|9626243|ref|NC_001416.1| 6 G N DP10;AF0.75534 gi|9626243|ref|NC_001416.1| 7 C N DP10;AF0.75535 gi|9626243|ref|NC_001416.1| 8 G N DP10;AF0.75536 gi|9626243|ref|NC_001416.1| 9 A N DP10;AF0.75537 gi|9626243|ref|NC_001416.1| 10 C N DP10;AF0.75538 539#### View MORE information with bcftools540 541 $ bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%FILTER[\t%TGT\t%DP]\n' sample1.fastq.gz-reference.vcf.gz | head542 gi|9626243|ref|NC_001416.1| 1 G N DP10;AF0.75 N/N 2543 gi|9626243|ref|NC_001416.1| 2 G N DP10;AF0.75 N/N 4544 gi|9626243|ref|NC_001416.1| 3 G N DP10;AF0.75 N/N 4545 gi|9626243|ref|NC_001416.1| 4 C N DP10;AF0.75 N/N 4546 gi|9626243|ref|NC_001416.1| 5 G N DP10;AF0.75 N/N 4547 gi|9626243|ref|NC_001416.1| 6 G N DP10;AF0.75 N/N 5548 gi|9626243|ref|NC_001416.1| 7 C N DP10;AF0.75 N/N 5549 gi|9626243|ref|NC_001416.1| 8 G N DP10;AF0.75 N/N 6550 gi|9626243|ref|NC_001416.1| 9 A N DP10;AF0.75 N/N 6551 gi|9626243|ref|NC_001416.1| 10 C N DP10;AF0.75 N/N 6552 553#### View the definitions for FILTER codes554 555For example, `DP10` indicates that the user set 10x as a threshold and the site has less than 10x coverage.556 557 $ zgrep ""##FILTER"" sample1.fastq.gz-reference.vcf.gz558 ##FILTER=<ID=DP10,Description=""Depth is less than 10, the user-set coverage threshold"">559 ##FILTER=<ID=RF0.75,Description=""Reference variant consensus is less than 0.75, the user-set threshold"">560 ##FILTER=<ID=AF0.75,Description=""Allele variant consensus is less than 0.75, the user-set threshold"">561 ##FILTER=<ID=isIndel,Description=""Indels are not used for analysis in Lyve-SET"">562 ##FILTER=<ID=masked,Description=""This site was masked using a bed file or other means"">563 ##FILTER=<ID=str10,Description=""Less than 10% or more than 90% of variant supporting reads on one strand"">564 ##FILTER=<ID=indelError,Description=""Likely artifact due to indel reads at this position"">565 566### Find constant sites567 568Lyve-SET focuses on variable sites, but what if you wanted to look at how many constant sites there are?569This is at least useful for when you need to tell BEAST or other similar programs what your constant sites are.570 571To get the full alignment, you can run the following `bcftools` command, 572followed by some parsing573 574```bash575bcftools query -f '%CHROM\t%POS\t%REF\t[%TGT\t]\n' --print-header out.pooled.vcf.gz > full.tsv.unrefined.tmp576 577# Get the first nucleotide of every nucleotide call to convert from diploid to haploid578perl -lane '579 BEGIN{580 # print the header581 $line=<>;582 chomp($line);583 print $line;584 $numFields=scalar(split(/\t/,$line));585 $lastIndex=$numFields-1;586}587for(@F[3..$lastIndex]){588 $_=substr($_,0,1);589}590print join(""\t"",@F);591' < full.tsv.unrefined.tmp > full.tsv592# Now you have full.tsv593 594# Parse full.tsv to get counts of constant sites595tail -n +2 full.tsv | perl -MData::Dumper -lane '596 # Print a dot every 100k lines597 print STDERR ""Looked at $. lines"" if($. % 100000 == 0);598 599 # Remove the contig, pos, ref fields600 splice(@F,0,3);601 # Pretend the first sample is the reference602 # and remove it from the list of samples with shift()603 $ref=shift(@F);604 # Skip if ""reference"" base is ambiguous605 next if($ref eq ""N"");606 # The site is constant until proven otherwise607 $is_constant=1;608 # Loop through all samples after the first609 for my $nt(@F){610 # If the samples nucleotide is not equal to the first samples,611 # then label as not constant and go to the next site612 if($nt ne $ref){613 $is_constant=0;614 last;615 }616 }617 $const{$ref} += $is_constant;618 END{619 print Dumper \%const;620 print STDERR ""\n"";621 }622'623```624 625This will give you output similar to626 627```text628$VAR1 = {629 'A' => 638550,630 'T' => 683152,631 'C' => 443617,632 '.' => 0,633 'G' => 412973634 };635```636 637Where you get counts for every nucleotide where it was constant at any given site in the set of samples.638You might be able to put it into BEAST with some syntax similar to (but not necessarily):639 640```xml641 <alignment dataType=""nucleotide"">642 <sequence idref=""taxon1"" value=""G"" count=""412973""/>643 <sequence idref=""taxon1"" value=""C"" count=""443617""/>644 <sequence idref=""taxon1"" value=""T"" count=""683152""/>645 <sequence idref=""taxon1"" value=""."" count=""0""/>646 <sequence idref=""taxon1"" value=""A"" count=""638550""/>647 <!-- Add more sequence elements for each taxon and its corresponding nucleotide with their respective counts -->648 </alignment>649```650 651And then, lastly, to get the full alignment including constant sites, you can transform `full.tsv` like so652 653```bash654matrixToAlignment.pl full.tsv > full.fasta655```656 657## A word on Grapetree and other minimum spanning tree software658 659It has been asked before if you can use Grapetree or something similar with Lyve-SET results.660First, yes you can but second, _why_?661A minimum spanning tree (MST) is useful for grouping similar profiles into the same circle.662Circles of profiles are linked by lines of a certain distance.663Then, circles of profiles are larger or smaller depending on how many members are in the circle.664This is a great way to visualize something like an outbreak.665However, it is not so great if every single sample has a different profile.666If you have different profiles, suddenly the MST becomes less informative and more chaotic.667 668If you still want to do this, then here are some example steps on how to run Grapetree on Lyve-SET results.669 670```bash671# make the profile spreadsheet from the alignment 672# by formatting the alignment into two-lines-per-entry fasta673seqtk seq -l 0 out.informative.fasta | \674 perl -lane '675 # Get the defline676 $sample=$_; 677 # Get the sequence678 $seq=<>; 679 # Remove the newlines680 chomp($sample, $seq);681 # Remove the > from the defline 682 $sample =~ s/>//; 683 # Transform the sequence into a set of sites in an array684 @seq=split(//, $seq); 685 # Print the profile 686 print join(""\t"", $sample, @seq);687 # Find out how many sites there are for the next step using STDERR 688 print STDERR ""# numSites: "".scalar(@seq);689 ' > profile.tmp.tsv690# numSites: 168691# numSites: 168692# numSites: 168693# numSites: 168694```695 696Now that we have a profile in `profile.tmp.tsv`, we still need a header using `168` sites.697 698```bash699# Generate a header of sites. ""0"" for the column of samples.700seq 0 168 | tr '\n' '\t' > profile.tsv701# Punctuate the header with a newline702echo >> profile.tsv703# Grab the data704cat profile.tmp.tsv >> profile.tsv705# Run grapetree however you want. Here is a very simple invocation.706grapetree --profile profile.tsv707(sample1:85,sample2:80,sample3:75,sample4:0);708```709 710## Other manual steps in Lyve-SET711 712### From a set of VCFs to finished results713 714 mergeVcf.sh -o msa/out.pooled.vcf.gz vcf/*.vcf.gz # get a pooled VCF715 cd msa716 pooledToMatrix.sh -o out.bcftoolsquery.tsv out.pooled.vcf.gz # Create a readable matrix717 filterMatrix.pl --noambiguities --noinvariant < out.bcftoolsquery.tsv > out.filteredbcftoolsquery.tsv # Filter out low-quality sites718 matrixToAlignment.pl < out.filteredbcftoolsquery.tsv > out.aln.fas # Create an alignment in standard fasta format719 set_processMsa.pl --numcpus 12 --auto --force out.aln.fas # Run the next steps in this mini-pipeline720 721### Manual steps in `set_processMsa.pl`722 723Hopefully all these commands make sense but please tell me if I need to expound.724 725 cd msa726 removeUninformativeSites.pl --gaps-allowed --ambiguities-allowed out.aln.fas > /tmp/variantSites.fasta727 pairwiseDistances.pl --numcpus 12 /tmp/variantSites.fasta | sort -k3,3n | tee pairwise.tsv | pairwiseTo2d.pl > pairwise.matrix.tsv && rm /tmp/variantSites.fasta728 set_indexCase.pl pairwise.tsv | sort -k2,2nr > eigen.tsv # Figure out the most ""connected"" genome which is the most likely index case729 launch_raxml.sh -n 12 informative.aln.fas informative # create a tree with the suffix 'informative'730 applyFstToTree.pl --numcpus 12 -t RAxML_bipartitions.informative -p pairwise.tsv --outprefix fst --outputType averages > fst.avg.tsv # look at the Fst for your tree (might result in an error for some trees, like polytomies)731 applyFstToTree.pl --numcpus 12 -t RAxML_bipartitions.informative -p pairwise.tsv --outprefix fst --outputType samples > fst.samples.tsv # instead of average Fst values per tree node, shows you each repetition732 733","Markdown"
734"Microbiology","lskatz/lyve-SET","docs/FAQ.md",".md","3025","40","Reference genomes735=================736 737Can I have multiple reference genomes?738--------------------------------------739No. However, you can add multiple assemblies to the asm directory. Assemblies don't have the same error profile as reads and so you might expect some skew in the ultimate phylogeny.740 741How do I choose a reference genome?742-----------------------------------743The best reference genome for an outbreak is something in-clade. You might even need to assemble the genome from your own reads before starting. An outgroup is not as related to your clade by definition and so it is probably not the best genome to use. The next best quality is a closed genome, or a genome with a high N50.744 745How to I include my reference genome in my analysis?746----------------------------------------------------747You might have noticed that your reference genome is not included in the final analysis.748If you want to include the reference genome in your final tree or SNP matrix, there are about two ways.749Either:750 751* Copy the assembly into the asm subfolder. There is a more automated way to do that with the `set_manage.pl` script, using `--add-assembly` as shown [here](EXAMPLES.md#prepare-the-project-directory)752* Find the original raw reads, e.g. Illumina fastq files, and place them into he reads folder. You can directly copy the interleaved files into that folder, or you can run `set_manage.pl --add-reads` on those interleaved files.753 754**NOTE** Keep the assembly in the `reference/` subfolder and continue to point `-ref` to the assembly in the `reference` subfolder.755 756High-quality-ness757=================758 759What are the different ways that SNPs in Lyve-SET have high confidence?760-----------------------------------------------------------------------761 762High quality for SNPs indicates that the resulting phylogeny will be high-fidelity. Although some SNPs are discarded that we are less sure about, the SNPs that we _are_ most sure about are retained, and the resulting phylogeny is the best inference.763 764User Lori Gladney created [this flowchart](../images/Lyve-SET_masking_mindmap_11-20-17.pdf) to help understand these points.765 766* **Detection of troublesome regions** such that they are not considered in hqSNP analysis. Currently in v1.0, only phage genes are detected; however other databases could be added in the future, and also I am open to other suggestions. Users can also specify a BED-formatted file to describe regions to mask.767* Only **unambiguous mapping** allowed768* Default **75% consensus** and **10x** coverage thresholds. These options can be changed when you launch Lyve-SET.769* Mechanisms to **remove clustered SNPs** -- not on by default however.770* **Maximum likelihood** phylogeny reconstruction. Ascertainment bias is also considered through RAxML v8.771* Both **forward and reverse reads** must support each SNP.772* Each read must have **95% identity** to the reference genome. In other words, if the read is 100bp, then only 5 differences in the read can be tolerated before the read is discarded.773","Markdown"
774"Microbiology","lskatz/lyve-SET","docs/VIZ.md",".md","1794","22","# Visualization775 776A large question is, after you are finished with a Lyve-SET run, how do you visualize the results? One of the advantages of Lyve-SET is its use of standard file formats. Therefore most files can be visualized in standard software. All output files are documented under [OUTPUT.md](OUTPUT.md).777 778## Tree779 780You can visualize the tree in many different tree drawing programs out there. Some of my personal favorites are [MEGA](http://www.megasoftware.net) and [Figtree](http://tree.bio.ed.ac.uk/software/figtree). When prompted, open `tree.dnd`. Some commercial software includes BioNumerics, CLC, and Geneious.781 782## SNP positions783 784SNP positions are encoded in the vcf files under the vcf directory. You can view them in any standard viewer such as [IGV](http://software.broadinstitute.org/software/igv). On the command line, although difficult to visualize, you can use `bcftools` which is redistributed in Lyve-SET.785 786To get started on IGV, first load the reference genome assembly. Second, load the bam and vcf files. You will immediately be able to browse the genome with these bam and vcf tracks. However for the advanced features, there is a learning curve (but it is worth it).787 788## Read alignments789 790Read alignments are encoded in bam files under the bam directory. These can be viewed with [IGV](http://software.broadinstitute.org/software/igv) and many commercial software packages such as CLC and Geneious. Additionally if you are command-line-inclined, you can view them with `samtools tview` which is redistributed with Lyve-SET.791 792## SNP distances793 794SNP distances might be the easiest to visualize. These can be viewed in LibreOffice or Microsoft Excel. Simply open `out.pairewiseMatrix.tsv` and turn on conditional formatting. This shows the distance heatmap.795","Markdown"
796"Microbiology","lskatz/lyve-SET","docs/INSTALL.md",".md","2697","89","Installation797============798 799## Containers800 801Here are the methods for installing with either Docker or Singularity.802 803 docker pull staphb/lyveset:1.1.4f804 805 singularity build lyveset.1.1.4f.sif docker://staphb/lyveset:1.1.4f806 807Requirements808------------809 810* **Perl, multithreaded**811 * BioPerl812* **BLAST+**813* **Edirect** (for downloading of test data)814* **GIT**, **SVN** (for installation and updating)815 816### Other requirements817 818Usually these packages are installed, but if you have a totally fresh system, you might want to consider installing the following819 820* zlib821* curses822* zip/unzip823* build-essential824 825Some Perl modules are needed in earlier versions of Lyve-SET including v1.1.4f.826 827* File::Slurp828* URI::Escape829 830Quickie Installation831--------------------832 8331. Run `make install` while you are in the Lyve-SET directory. This step probably takes 10-20 minutes.8342. Update the path to include the scripts subdirectory. You can do this yourself if you are comfortable or run `make env`.8353. Update your current session's path: `source ~/.bashrc`836 837In-depth Installation838---------------------839 840Download the latest stable version from https://github.com/lskatz/lyve-SET/releases. You can also roll the dice by getting the cutting edge version with git.841 842 tar -zxvf lyve-SET-1.1.4f.tar.gz843 cd lyve-SET-1.1.4f.tar.gz844 make install845 make env846 847Other Installation Options848------------849* `make install`850* `make install-optional` to install optional prerequisite software including bwa and snap851* `make env` - update `PATH` and `PERL5LIB` in the `~/.bashrc` file.852* `make check` - check and see if you have all the prerequisites853* `make test` - run a test phage dataset provided by CFSAN854* `make help` - for other `make` options855* `make clean` - clean up an old installation in preparation for a new installation856 857### Fine details858 859* `make install-*` - Many other installation options are available including but not limited to:860 * `make install-smalt`861 * `make install-CGP`862 * `make install-samtools`863* `make clean-*` - Every `make install` command comes with a `make clean` command, e.g.:864 * `make clean-CGP`865 866Upgrading867---------868 869### By stable releases870Unfortunately the best way to get the next stable release is to download the full version like usual, followed by `make install`. If successful, then delete the directory containing the older version.871 872 cd ~/tmp873 wget https://github.com/lskatz/lyve-SET/archive/v1.1.4f.tar.gz874 tar zxvf release.tar.gz875 cd Lyve-SET876 make install # takes 10-20 minutes to download packages on broadband; install877 cd ~/bin878 rm -r Lyve-SET && mv ~/tmp/Lyve-SET .879 880### By `git`881 git pull -u origin master882 make clean883 make install884","Markdown"
885"Microbiology","lskatz/lyve-SET","docs/TROUBLESHOOTING.md",".md","5827","75","# Troubleshooting886 887Things do not always go smoothly. Here is a list of things to try out if things go wrong.888 889## Interleaved reads890 891Why are both R1 and R2 showing up in my results?892 893Lyve-SET requires _interleaved_ reads (or commonly called ""shuffled""). This means that there is one fastq file per genome and that their format is such that there is one forward read, followed by its reverse reads, followed by the next forward read, etc.894To interleave one set of reads, you can use the script `run_assembly_shuffleReads.pl`. To shuffle many pairs of reads, you can use `shuffleSplitReads.pl`. Both have usage statements if you run them with `--help`.895 896## Smalt says that there is an invalid fastq/fasta format897 898This can be the result of a few different things. 899 900### Reference assembly integrity901 902First, check to see if your reference genome is intact.903 904* Manually inspect the fasta file with `less`, e.g., `less reference/reference.fasta`.905* Another way to manually inspect the fasta file is with `grep`, e.g., `grep -A 1 "">"" reference/reference.fasta`.906* Get some summary metrics with `run_assembly_metrics.pl`, e.g., `run_assembly_metrics.pl reference/reference.fasta | column -t`.907 908### Fastq integrity909 910You can also see if your fastq files are intact which is a little more tricky.911 912* Are your files in the interleaved format? See [Interleaved reads](#interleaved-reads).913* Are your files gzipped? Use the Linux command `file` to see that it says 'gzip' in the description of your files. For example, `file reads/*.fastq.gz`. It is not necessary to have gzipped files, but it _is_ necessary for the extension to match the type. For example you do not want to have compressed files that end in `.fastq.gz` or uncompressed files that end in `.fastq`.914* Are the actual files intact? 915 * You can use the [lskScript](https://github.com/lskatz/lskScripts) tool `validateFastq.pl` to find common mistakes. For example: `zcat reads/genome.fastq.gz | validateFastq.pl --pe --min-length 1 --verbose`.916 * Is the interleaved file line count the number of lines of both R1 and R2? Example command: `zcat R1.fastq.gz | wc -l`.917 * Simply reshuffle: `shuffleSplitReads.pl -o shuffled --numcpus 1 split/*.fastq.gz`918 919## BCFtools merge920 921### cannot read fastq header; error with bcftools merge922 923This is likely due to [Fastq integrity](#fastq-integrity)924 925## The tree was not created926 927### Too few taxa928 929RAxML requires at least four genomes to create a tree. If you have few genomes in your analysis, perhaps the `out.pairwise.tsv` or `out.pairwiseMatrix.tsv` file would make more sense in terms of deliverable results.930 931### Nonunique names932 933* RAxML only considers the **first 30 characters** of your genome name. In the offhand chance that you have two genomes with filenames that vary on the 31st or later characters, then consider renaming your files appropriately.934* Did you include **R1 and R2 and/or also the shuffled reads**? Lyve-SET requires shuffled reads and so if you haven't done this, please follow the steps under the section [Interleaved reads](#interleaved-reads). You should not be using both R1 and R2 in your input directly. If they are already shuffled, delete all instances of R1 and R2 files in the project directory. However, *do not* delete the shuffled file. *Do not* delete your original data (just a general tip). Delete R1/R2 files and their derivative files in the reads subdirectory, the tmp subdirectory, the bam subdirectory, and the vcf subdirectory. Delete all files in the msa directory, because they will have to be recreated. Rerun Lyve-SET. This is a huge issue because935 1. R1 and R2 SNPs will not be as accurate as the paired-end shuffled entry, because their read mapping is not as accurate936 2. These less accurate entries will enter noise into the tree inferrence.937 938### Not enough variable sites939 940#### Best case scenario941 942This isn't a huge problem but it's not immediately obvious that this is the case always. If all of your genomes have the exact same SNP profile, then the multiple sequence alignment will have zero variable sites, and RAxML will fail. This is a result in itself because it says that all your genomes have the same genomic profile.943 944How do you detect this? You can view `out.snpmatrix.tsv` to see if there are many sites without majority ambiguous allele calls. You can see if `out.snpmatrix.tsv` has many sites. Therefore you can do a simple calculation to show that you have many high-quality sites but few or no variable sites.945 946How do you avoid this situation? If you have a good outgroup, then there will be variable sites, and it won't matter if your clade of interest has no variation between them. RAxML can build a tree when there are a few guaranteed variable sites introduced by the outgroup. This outgroup by definition is a genome that is phylogenetically related but is not the same profile as your clade of interest.947 948#### Worst case scenario949 950What if they are *supposed* to have different SNP profiles?951 9521. Check to see if the **number of masked sites is huge** in the multiple sequence alignment. `launch_set.log` in the log directory should inform you of huge percentages of masked sites.9532. Check to see if there are weird read mappings. See the [visualization guide](VIZ.md) for more details. Or some basic commands954 * `samtools depth bam/somegenome.bam | less`955 * `samtools flagstat bam/somegenome.bam`956 * `bam stats --basic --in bam/somegenome.bam`9573. Check to see if there are a ton of **failed sites** in the VCF files. This is more complicated but it will give you a trend. For example, do you see a lot of low coverage sites? Low consensus sites?958 * Example command to view the number of sites passing, or number of sites failed per reason `bcftools view --no-header somegenome.vcf.gz | cut -f 7 | sort | uniq -c`959","Markdown"
960"Microbiology","lskatz/lyve-SET","docs/EXAMPLES.md",".md","6596","148","Examples961========962 963Run a test dataset964------------------965 966The script `set_test.pl` will run an actual test on a given dataset. It uses `set_downloadTestData.pl` to get any bacterial genomes and then runs `launch_set.pl`. However, the lambda dataset is small enough to fit on GIT and does not need to be downloaded.967 968 Runs a test dataset with Lyve-SET969 Usage: set_test.pl dataset [project]970 dataset names could be one of the following:971 escherichia_coli, lambda, listeria_monocytogenes, salmonella_enterica_agona972 NOTE: project will be the name of the dataset, if it is not given973 974 --numcpus 1 How many cpus you want to use975 --do-nothing To print the commands but do not run system calls976 977`$ set_test.pl lambda # will run the entire lambda phage dataset and produce meaningful results in ./lambda/msa/`978 979 980Prepare the project directory981-----------------------------982 983The script `set_manage.pl` sets up the project directory and adds reads, and you should use the following syntax:984 985 $ set_manage.pl --create setTest986 987Depending on your knowledge of Linux, you might choose to set up the rest of the project using `set_manage.pl` or using symlinks. This is the `set_manage.pl` way:988 989 $ set_manage.pl setTest --add-reads file1.fastq.gz990 $ set_manage.pl setTest --add-reads file2.fastq.gz991 $ set_manage.pl setTest --add-reads file3.fastq.gz992 $ set_manage.pl setTest --add-reads file4.fastq.gz993 $ set_manage.pl setTest --add-reads file5.fastq.gz994 $ set_manage.pl setTest --add-assembly file1.fasta995 $ set_manage.pl setTest --add-assembly file2.fasta996 $ set_manage.pl setTest --change-reference file3.fasta997 998This is the symlink way:999 1000 $ cd setTest/reads1001 $ ln -sv path/to/reads/*.fastq.gz . # symlink reads1002 $ cd ../asm1003 $ ln -sv path/to/assemblies/*.fasta . # symlink the assemblies1004 $ cd ../reference1005 # copy the assembly in case you need to alter it (e.g., remove small contigs or edit deflines)1006 $ cp -v path/to/assemblies/reference.fasta .1007 1008Run Lyve-SET with as few options as possible1009 1010 $ launch_set.pl setTest1011 1012More complex1013 1014 $ launch_set.pl setTest --queue all.q --numnodes 20 --numcpus 12 --noclean --notrees1015 1016If you specified notrees, then you can edit the multiple sequence alignment before analyzing it. See the next section on examples on how/why you would edit the alignment.1017 1018 $ cd setTest/msa1019 $ gedit out.aln.fas # alter the deflines or whatever you want before moving on1020 # => out.aln.fas is here1021 $ set_process_msa.pl --auto --numcpus 121022 # Optionally, qsub this script instead because it could be cpu-intensive1023 $ qsub -pe smp 12 -cwd -V -o trees.log -j y -N msaLyveSET -S $(which perl) $(which set_process_msa.pl) --auto --numcpus 121024 1025If you make a mistake and need to redo something:1026 1027 # remove all intermediate files1028 $ rm setProj/bam/genome*1029 $ rm setProj/vcf/genome* setProj/vcf/unfiltered/genome*1030 # OR, remove a genome entirely1031 $ set_manage.pl setProj --remove-reads genome.fastq.gz1032 $ set_manage.pl setProj --remove-assembly genome.fasta1033 # remove the last multiple sequence alignment files1034 $ rm -r setProj/msa/*1035 # or save the MSA results for another time1036 $ mv setProj/msa setProj/msa.oldresults && mkdir setProj/msa1037 # redo the analysis (all untouched bams and vcfs will not be redone)1038 $ launch_set.pl setProj1039 1040Why would you want to edit the out.aln.fas file? Or what kinds of things can you observe here before making a tree?1041 1042 # Alter the identifiers of your genomes, so that they look nice in the phylogeny(ies)1043 $ sed -i.bak 's/\.fastq\.gz.*//' out.aln.fas1044 # OR maybe to just remove the reference genome name:1045 $ sed -i.bak 's/\-reference.*//' out.aln.fas1046 1047 # Be sure that the taxa names are unique still.1048 # If there is any output, then you have duplicated names which need to be fixed.1049 $ grep "">"" out.aln.fas | sort | uniq -d # nothing should show up1050 1051 # After the taxon names are edited nicely,1052 # back up the out.aln.fas file before any entries are removed1053 $ cp -v out.aln.fas out.aln.fas.bak1054 1055 # => All extensions are removed in taxon names; a backup of the file was named out.aln.fas.bak1056 # find the genomes with the most number of Ns (ie masked SNP calls)1057 $ perl -lane 'chomp; if(/^>/){s/>//;$id=$_;}else{$num=(s/(N)/$1/gi); print ""$id\t$num"";}' < out.aln.fas|sort -k2,2n|column -t1058 # => consider removing any genome with too many masked bases by manually editing the file.1059 1060 # Run SET on your new set of genomes out.aln.fas.1061 $ set_process_msa.pl out.aln.fas --auto --numcpus 121062 1063 # Create a new subset as needed, but you should always read from your master record and not the informative.aln.fas file1064 # which is now out.aln.fas.bak1065 $ cp -v out.aln.fas.bak out.aln.fas1066 $ ... # more editing here1067 $ set_process_msa.pl out.aln.fas --auto --numcpus 121068 1069 1070Specific ways to regenerate files1071---------------------------------1072 1073Lyve-SET is very modular and so there are specific scripts to regenerate files. Since you might edit the msa files, you might want to know how to recover in case you make a mistake.1074 1075Need to remake out.aln.fas?1076 1077 $ mvcfToAlignment.pl out.pooled.vcf.gz --bcfOutput bcfquery.out > out.aln.fas1078 1079Need to remake out.pooled.vcf.gz? Use bcftools.1080 1081 $ bcftools merge vcf/unfiltered/*.vcf.gz -O -z > msa/out.pooled.vcf.gz1082 $ tabix -f msa/out.pooled.vcf.gz # Always index your compressed vcf files1083 1084Need to remake all the other msa files after you recreated out.aln.fas?1085 1086 $ set_processMsa.pl --auto out.aln.fas --numcpus 12 --force1087 1088Need to make a backup of your current results? Just rename the directory but make sure that the msa directory still exists.1089 1090 $ mv msa/ msa.bak/ && mkdir msa1091 1092Special cases1093-------------1094 1095Q: What if you have a non-Illumina file that you would like to use?1096 1097A: You can still map your reads separately into a sorted bam file and then run SET after that bam file has been created. For example, if you have a 454 or Ion Torrent file, you can map it using bwa-sw with the following steps. Keep the file naming system intact.1098 1099 bwa index -a bwtsw reference/2010EL-1786.fasta1100 bwa bwasw -t 12 reference/2010EL-1786.fasta reads/example.fastq.gz > tmp/example.fastq.gz.sam1101 samtools view -bS tmp/example.fastq.gz.sam > tmp/example.fastq.gz.bam1102 samtools sort tmp/example.fastq.gz.bam bam/example.fastq.gz-2010EL-1786.sorted # produces the bam extension1103 samtools index bam/example.fastq.gz-2010EL-1786.sorted.bam1104 rm -v tmp/example*1105 1106 launch_set.pl ...1107","Markdown"
1108"Microbiology","lskatz/lyve-SET","scripts/launch_phyml.sh",".sh","906","38","#!/bin/sh1109#$ -S /bin/sh1110#$ -pe smp 11111#$ -cwd1112#$ -V1113#$ -o launch_phyml.sh.out -j y1114 1115# makes a tree out of an aln1116 1117script=`basename $0`;1118aln=$11119if [ ""$aln"" = """" ]; then1120 echo Usage: `basename $0` aln.phy1121 exit 1;1122fi1123 1124# find which phyml to use1125phyml=`(which phyml-mpi || which phyml || which PhyML || which PhyML-3.1_linux64 || which phyml_linux_64) 2>/dev/null`;1126if [ $? -gt 0 ]; then echo ""Could not find phyml""; exit 1; fi;1127echo ""$script: Found phyml at $phyml""1128 1129# get the extension1130b=`basename $aln`;1131suffix=""${b##*.}"";1132if [ ""$suffix"" != ""phy"" ]; then1133 echo ""Converting to phylip format because I did not see a phy extension"";1134 convertAlignment.pl -f phylip -i $aln -o $aln.phy;1135 if [ $? -gt 0 ]; then exit 1; fi;1136 aln=""$aln.phy"";1137fi;1138 1139 1140yes | \1141 $phyml -i $aln -b -4 -m GTR -s BEST --quiet;1142if [ $? -gt 0 ]; then echo ""$script: ERROR in phyml"" exit 1; fi;1143echo ""$script: Finished without error!"";1144 1145","Shell"
1146"Microbiology","lskatz/lyve-SET","scripts/testInstallation.sh",".sh","112","9","#!/bin/bash1147 1148set -e1149 1150set_test.pl --numcpus 2 lambda lambda -- --noqsub --nodiagnose1151 1152echo ""Everything passed!""1153 1154","Shell"
1155"Microbiology","lskatz/lyve-SET","scripts/pooledToMatrix.sh",".sh","1272","69","#!/bin/bash1156# Runs bcftools query to make a SNP matrix1157 1158script=$(basename $0);1159 1160usage () {1161 echo ""$script: generates a SNP matrix using BCFtools and a pooled VCF file""1162 echo ""USAGE: $script -o bcfmatrix.tsv pooled.vcf.gz""1163 return 0;1164}1165 1166logmsg () {1167 echo ""$script: $@"";1168 return 0;1169}1170 1171while getopts ""ho:"" o; do1172 case ""${o}"" in1173 h)1174 usage1175 exit 11176 ;;1177 o)1178 OUT=""$OPTARG""1179 ;;1180 *)1181 echo ""ERROR: I do not understand option $OPTARG""1182 ;;1183 esac1184done1185shift $(($OPTIND-1)) # remove the flag arguments from ARGV1186 1187IN=$11188 1189if [ ""$OUT"" == """" ] || [ ""$IN"" == """" ]; then1190 usage1191 exit 1;1192fi1193 1194# Create the basic matrix1195t='\t'1196command=""bcftools query -i '%TYPE=\""snp\""' -f '%CHROM$t%POS$t%REF$t[%TGT$t]\\n' --print-header $IN > $OUT.unrefined.tmp""1197logmsg $command;1198eval $command1199if [ $? -gt 0 ]; then1200 echo -e ""ERROR with bcftools:\n $command"";