dataesr/scientific-paragraphs-categorization huggingface dataset:
There are several large consortium genome sequencing projects producing genomic data at an ever increasing rate (e. g. ,  21, 22) , providing raw data for SNP discovery. Even when a complete, high-coverage genome assembly for a species of interest is not available, generating low-coverage genome data and aligning to a related species can provide relatively quick and inexpensive access to thousands of SNPs in a target species. The methods in this chapter assume availability of a Unix or Linux computing system with sufficient capacity (memory and storage) to handle Gigabyte data sets. The primary software packages are all freely available, though some additional commercial software is recommended (but not necessary) for some analyses.
Materials The methods in this chapter are based on the versions of software listed in Table  1. There is no one public software package that does all of the steps described below, so we have used multiple freely-available programs for data processing. We also describe several scripts that have been written using different programs (e. g. , Python, R, Perl), so the programs to run these scripts are necessary. Software can change, so it is important to use the specified versions, or test newer versions to make sure that the commands are performing the same functions and options remain as stated. The best way to do this (and familiarize yourself with the software) is to read the help documents or guidelines for each software package and the specified commands. The current versions of the scripts described below are available through GitHub at https://github.
If there are multiple biosamples, it is important to determine which SRA files are from each individual if the goal is to use either pooled sample data or just a single individual for SNP discovery. On the BioProject description web page, the Assembly link leads to the project summary. All of the assembly files can be obtained via ftp through the NCBI site: ftp://ftp. ncbi. nlm. nih. gov/genomes/. After connecting as a guest, you can navigate to the species directory and view the README file that describes each of the directories and files within. For our purposes, the important files are in the directory CHR_Un, which contains the "Chromosome directories", or "concatenated sequence data for scaffolds that have been assembled from individual GenBank records", in a variety of formats. The files ending in ".
fa. gz" are contigs in compressed FASTA format, and the ones ending in ". mfa. gz" are repeat-masked contigs in compressed FASTA format. These can be downloaded by a variety of methods. Probably the most useful in a Linux environment is to use "wget", for example: File transfer may take several hours, depending on the size of the file and speed of connection. It is usually not advised to map reads to a repeat masked genome, as reads neighboring the masked repeats may not map, thereby reducing coverage around the repeats. However, if you have a reference genome with reasonably high coverage, the loss will be minimal relative to the amount of high-coverage genome data, and using the repeat-masked reference sequence will reduce the need for filtering out SNPs in repeats later. The next step is to identify and download the short-read archive data.
The files can be selected by checking the boxes of those you want and using the "Download" button. This generates a web page with the download buttons for individual files, or you can select files and click the button for "download SRA URLs" to download a text file containing all of the URL links to facilitate command-line transfer of the files using the "wget" command in a Linux environment: Example SRA URLs file: download_sra_urls. txt (for the Orcinus orca project): ftp://ftp-trace. ncbi. nlm. nih. gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR574/SRR574967/SRR574967. sra ftp://ftp-trace. ncbi. nlm. nih. gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR574/SRR574968/SRR574968. sra ftp://ftp-trace. ncbi. nlm. nih. gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR574/SRR574969/SRR574969. sra (etc. ) # Transfer each of the files using "wget" and the URLs in the # download_sra_urls.
This requires the fastq-dump function from the SRAtoolkit software. Files size will expand substantially, typically by 4-5x the original SRA file size. # Convert each SRA files into FASTQ files fastq-dump SRR574967. sra --split-3 The "--split-3" option indicates that the file contains paired-end reads and will be split into two fastq files, one for read 1 and a second for read 2 (e. g. , in this case, SRR574967_1. fastq, SRR574967_2. fastq). If the file only contains single-read data, only a single file will be generated, so it is recommended that this option is always used to ensure data are correctly parsed into files. The FASTQ file may be too large to be processed when only one SRA file is available on the SRA portal for a high coverage genome. Fastq-splitter can be used to divide the large file in several smaller files that are easier to handle.
After downloading the program, you need to make sure that the perl script is set as an executable. # Command to check if the perl script is set to be executable chmod +x fastq-splitter. pl # Split the file in 6 parts for both the forward and reverse reads. /fastq-splitter. pl --n-parts 6 --check SRAfileF. fastq. /fastq-splitter. pl --n-parts 6 --check SRAfileR. fastq # The path to the perl script (fastq-splitter. pl) should be given. # --n-parts indicates into how many files the large # SRA file should be divided (6 in this example). # --check does some verification of FASTQ format correctness.
Alignment of reads to reference contigs: Alignment of reads to reference contigs is fairly straightforward, but requires multiple steps that can be confusing. There are multiple alignment programs, but one of the most frequently used and easiest to use is the Burrows-Wheeler Aligner (BWA; 26). When there are multiple smaller FASTQ files for a single individual, they can be concatenated into a single file for each read direction. This is a Linux function: # Concatenate forward read files. cat filename1_read1. fastq filename2_read1. fastq > forward_reads. fastq However, when the files are large (e. g. , >10GB), running separate alignments of each FASTQ file (or pair of files for paired-end data) will improve speed by allowing parallel processes run on separate cores (see Note 1). The resulting alignments can be combined after conversion to the smaller BAM files (see below).
Although there are many options that can be altered to vary alignment parameters, the updated aligner "BWA-MEM" algorithm  (27)  is implemented in BWA and has been optimized for aligning short reads to a large genome such as human without the need to specify alignment parameters, and performs as well as or better than many other aligners with short read (100bp) and longer read-length data  (27). Mitochondrial DNA can make up a significant portion of sequence reads, so removal of reads that map to the mitogenome can reduce file sizes and speed up subsequent analyses. A reference mitogenome FASTA file can typically be downloaded from GenBank for the target species or a closely related species, or the mitogenome sequence can often be assembled de novo from the NGS reads (e. g. ,  28). Map the reads to the mitogenome with BWA, convert the SAM file to a BAM file and sort with SAMtools.
samtools index individual1_q30. bam The SNP discovery pipeline is managed through a shell script (see Appendix 1) that performs several functions. The shell script is launched by executing the script in a directory that also must contain the reference FASTA file (previously indexed with samtools faidx), the alignment BAM file (sorted) and the SNP discovery Python script "Script2_generate_genotype_blocks. py" (see Appendix 2). The shell script takes as input the alignment BAM file (sorted), reference FASTA file, and the total desired bases of the flanking sequence of the SNP (by default 300 bp). The shell script takes objects in order, assigning them to aliases that are used to fill in the appropriate filenames into a series of commands. Execution of the shell script first runs samtools mpileup to compile the read bases, quality and coverage information at each site in the alignment, then calls SNPs from the pileup using BCFtools, using filters to exclude loci with depth of coverage below and above specified levels (set within the shell script; see Appendix 1).
The custom Python script "Script2_generate_genotype_blocks. py" further filters the SNPs based on user-defined quality cut-off and flanking sequence length parameters, and outputs two results files (see Note 2). # 'output' represents the file name that will be assigned to the # output files. The output files from this script are a FASTA file containing the designated target sequences surrounding each SNP, and a genotype VCF file. The VCF file format is described in detail elsewhere (https://githum. com/samtools/hts-specs). This file contains information about each SNP, starting with the contig ID ("chrom"), the SNP position, the reference and alternate alleles, a quality score, and a list of additional "info" about the SNP, including the depth of coverage (DP), various quality measures, the count forward and reverse reads for each allele used in variant calling (DP4), and the consensus quality (FQ).
Sequence data processing This section assumes that single-or paired-end read data have been de-multiplexed to separate reads by the index sequences of each sample in the multiplex pool. SRA files downloaded from GenBank have presumably been processed to remove adapters and filtered for quality. For data generated specifically for SNP discovery, these steps need to be done prior to mapping reads to the reference genome. For each FASTQ file, ADAPTER-REMOVAL (or another program, e. g. , Trimmomatic; 32) can be used to remove adapter sequences from the sequencing reads and remove sequence reads that are ≤30 bp or another length threshold following trimming. Reads in FASTQ format must then be mapped to the reference genome of the studied species or a related species, then converted to a sorted BAM file following the same steps outlined above.
The reference can be repeat hard-masked, or repeats can be masked in the BAM file following mapping. To hard mask the reference FASTA file to which reads are being mapped to, we use RepeatMasker. Repeat libraries must also be accessed from the RepBase database which is available through the GIRI Web site (http://www. girinst. org/repbase/index. html). Note that this requires opening an account, however, there is currently no charge for this for academic non-profit users. The output files include a masked copy of the input file, which contains the query sequences, with identified repeats and low complexity sequences masked and replaced with 'N's. A map file is also generated, in which coordinates in terms of contig, start and end position are given, which can then be used as a. bed file (see below).
-Filter 7: (Optional) Remove SNPs with excess heterozygosity. This step can only be conducted meaningfully when there are sufficient samples to infer a significant deviation from the expectations of Hardy-Weinberg equilibrium (HWE), and is useful in identifying putative SNPs that are the result of gene duplications, gene families, and repetitive regions. -Filter 8: (Optional) Compare SNP locus sequences to reference databases (e. g. , GenBank) to filter out loci that might match to non-target species, duplicated loci, gene families, and repeat regions. The commands to run these 8 filtering steps are detailed below. Three output files are produced, including data. depthSample and data. depthGlobal. In data. depthSample each line represents one sample (i. e. , one individual) and column 1 is the number of sites with a depth of 0×, column 2 the number of sites with a depth of 1×, etc.
qc. arg, bam. qc. depthGlobal, bam. qc. depthSample and bam. qc. qs. These files can be plotted using the R script Script10_plotQC. R (see Appendix 12). The script creates a pdf file of the plotted output, and the file "bam. qc. info" with the q-score and global depth data. As an example, we consider here that the mean coverage across all samples is 5× and that regions with coverage >10× are potential unmasked repeated regions or mapping artifacts. Regions of poor mapping quality (Q<30) and excessive coverage (>10×) are detected using the CALLABLELOCI tool in GATK. Running CALLABLELOCI involves some data formatting. # Create a dictionary file on the initial reference using Picard tools. # -R = The reference genome against which the sequence data was # mapped. The GATK requires an index file (samtools faidx) and # a dictionary file accompanying the reference.
The SNP at position 601678 on scaffold JH472447 in the. mafs file is in a region of poor mapping quality (JH472447, positions 601658-601711). The R script is used to filter this SNP out. All the other SNPs in these example lines are not in a region of poor mapping in the example lines of the "poor_mapping_quality. txt" filter given above and will be kept. # Source the "remove_poor_mapping_quality" script (Appendix 5) # file using "Rscript" or execute individual lines from the # script in the R environment. Rscript Script4a_Filter1_Remove_poor_mapping_quality. R # This script requires the input files "genolike_data. mafs" # and "poor_mapping_quality. txt". # The output file is "SNPs_goodquality. txt" A similar script is used the remove the SNPs in regions of excessive coverage. The SNPs in those regions are likely to be in a repeat region, NUMT or in a region with some other mapping artifact.
Notes 1. If there are multiple sample files, running them sequentially through the same process can be facilitated by using shell scripts (e. g. , Script1_SNP_call_filter. sh, Appendix 1). A shell script simply a text document containing the commands that would be entered sequentially. For running the same commands on a set of sample files, the same command can be entered for each sample, with only the filenames (input and output) changed. The script is executed by typing '. /' before the script name in the directory where your script and input files reside. Permissions may need to be set appropriately to allow execution of shell scripts (see "chmod" in Linux commands, see Table  2 ). 2. Within the shell script text file, the parameters for minimum and maximum depth of coverage are set by the user, as determined for the average depth of coverage of the BAM alignment file.
The script also calls the python script, which must be in the same directory, and the path to the python program may also have to be specified. 3. The filenames provided in these steps are used later in R-scripts, so changing the names will require changes in the R-scripts as well. This process is slow, and may take days to complete. It does not typically use >1 thread, though it may increase thread use as it adds sample files. A newly published wrapper for the program ANGSD (40) may make some of these steps easier or improve visualization, but have not yet been tested by us. 4. If the mapping quality filter is applied and results in an empty output file (i. e. , SNPs with mapping quality <30 have all been previously filtered out), then the subsequent filter will produce an error. This can be circumvented by copying the first line of the "excessive_coverage.
txt" file and pasting it into the empty "poor_mapping_quality. txt" file prior to running Script4a. This will simply remove the one of the excessive coverage SNPs and allow the script to proceed. 5. Calculating mean depth of coverage from the file data. depthGlobal: Open the file in Excel, insert the depth of coverage from 0 to 100 in row 1 (above the depth data). Use a function to calculate the mean depth: =SUMPRODUCT(B1:CW1*B2:CW2)/ SUM(B2:CW2). B1:CW1 contain the different depth of coverage from 0 to 100 and B2:CX2 the number of sites for each depth of coverage. Tables: top Typing "top" shows the current processes on the system, including how the percent CPU (i. e. , 100% = 1 CPU) and percent of the total memory being used by each process. Exit with the keyboard combination "control" and "c" df -h shows current system status, including size, amount used, amount available, % used.
0. 6 http://www. repeatmasker. org/RMDownload. html (optional) SAMtools 1. 2 (45) SEQTK 1. 1-r92 https://github. com/lh3/seqtk SRAtoolkit 2. 5. 7 http://www. ncbi. nlm. nih. gov/Traces/sra/sra. cgi?view=soft ware/.
as a powerful approach, but the lack of sufficient genomic DNA material due to cells' unculturability prevents traditional shotgun sequencing and genome assembly. Currently, there still is a big knowledge gap in protist genomics caused by the reluctance of many species to grow in culture. Indeed, only 15% of the completed or ongoing projects in the Genomes OnLine Database (GOLD; https://gold. jgi. doe. gov)  15  concern protistan taxa, and most of them come from cultured phototrophic  [16] [17] [18]  or parasitic species  19  , resulting in a biased view of the full eukaryotic diversity  20. In this context, a natural option is single cell genomics (SCG), which produce Single Amplified Genomes (SAGs) that can later be sequenced. This technology was initially used to produce genomic information from single prokaryotic cells collected from the ocean  21.
As expected, the obtained topography (U-Matrix) representing the structure of the tetranucleotide frequency dataset formed two large clusters that coincided with MAST-4A and MAST-4E contigs (Fig. 3 ). All SAGs within the same lineage were found in the same region of the map, revealing the same tetranucleotide frequency profile. In addition, we observed two small clusters that contained particular genomic signatures: the subcluster "a" of mitochondrial origin and the subcluster "b" of putative prey origin. A detailed analysis of the contigs within these subclusters identified four subgroups (Fig. 3 ). The first (a1) included 17 contigs from both lineages related to the mitochondrion of Cafeteria roenbergensis (mean similarity of 94%); probably these belonged to MAST-4 mitochondrion. The three other subgroups derived from putative prey DNA.
9 Mbp  26. Compared with published genomes of picosized protists and other heterotrophic flagellates, the number of genes predicted in MAST-4 genomes is relatively high (Table  1 ), only comparable with the parasitic Phytophthora spp, and perhaps some could derive from foreign DNA. On the other hand, the mean gene size (1,657 and 1,723 bp in each lineage) was similar to the values found in other protists genomes (Table  1 ). Compared to other protists, MAST-4 has compact genomes with few but long introns (Table  1 ). Finally, a first gene annotation of the two co-assemblies were performed using BLASTp against the KEGG Orthology (KO) database 54 revealing a total of 2,733 and 2,210 good KO hits for MAST-4A and MAST-4E, respectively (Table  1 ). Similarly, predicted proteins were also searched against the eukaryotic orthologous groups (KOG) using rpsBLAST and the Conserved Domains and Protein Classification database of NCBI (NCBI-CDD)  55  and a total of 2,115 and 1,878 KOGs were assigned among the MAST-4A and MAST-4E co-assemblies (Table  1 ).
Sample associated environmental metadata are reported in Supplementary Table  S1  and more details can be found in PANGAEA  59. SAG sequencing, assembly, quality control and completeness assessment. After purification of the MDA products and generation of 101 bp paired-end libraries, each SAG was sequenced in a 1/8 th Illumina HiSeq lane at the Oregon Health & Science University (US) or the National Sequencing Center of Genoscope (France) (Supplementary Table  S2  and Supplementary Fig. S1a ). Reads were assembled or co-assembled using SPAdes 3. 1 or 3. 6  60. In all assemblies, contigs shorter than 1 kbp were discarded. Quality profiles and basic statistics (genome size, number of contigs, N50, GC content) of single SAG assemblies and co-assemblies were generated with Quast 51. Estimations of genome recovery were done with CEGMA 61 (Core Eukaryotic Genes Mapping Approach; Supplementary Fig.
S1d ). In order to assess if genome completeness in each SAG depended on sequencing effort, reads from the two largest SAG assemblies within each clade were randomly subsampled into 5 different sequencing depths with the seqtk toolkit (https://github. com/lh3/seqtk). Five independent replicates were generated for each sequencing depth using different random number generator seeds. For each pool of subsampled reads, new assemblies and genomes recoveries (contigs > 1 kbp) were generated as described above. SAG comparisons based on nucleotide identity and tetranucleotide frequency. Nucleotide identity between SAGs was estimated by a pairwise BLAST analysis 62 between full-length contigs of all SAGs within each clade, with a minimum similarity of 70% and a maximal e-value of 10 -5. Tetranucleotide frequencies in each individual SAG were calculated using a 1 bp sliding window in both DNA strands in contigs between 2.
5 and 5 kbp in size with a custom Perl script  43  and clustered using ESOM  42  (Emergent Self-organizing Maps; Supplementary Fig. S1b ). Raw data were normalized using robust estimates of mean and variance ("Robust ZT" option) and trained according to Dick and colleagues  43  with the k-Batch algorithm and Euclidean grid distance. Sub-clusters of interest were isolated to identify the corresponding contigs by a BLASTn analysis against the NCBI-nt and NCBI-RefSeq (including organelles genomes) databases  63. Blast hits (similarity > 80%, e-value < 10 -5 ) were taxonomically assigned. Genome analysis using fragment recruitment tools. Original reads were mapped back to their corresponding co-assembly using bowtie2 with default parameters 64 (Supplementary Fig. S1c ). The reads alignments (BAM file) obtained were processed using samtools  65  , BEDTools 66 , QualiMap2 67 and custom perl scripts.
Then a comparison of each individual SAG assembly against the co-assembly as a reference for both MASTs was performed using Quast  51. Analysis of Core Eukaryotic Genes (CEGs) and of the rDNA operon. A subset of 248 universal CEGs within each SAG and the two final co-assemblies were identified with CEGMA  61. For each detected CEG, amino acids sequences were aligned using Clustal-Omega  68. These alignments were then used to calculate distance matrices based on percent identities for each sequence pair. We searched for contigs containing the 18S rDNA sequence in all individual MAST-4 SAGs and the co-assembly. The complete rDNA operon sequences were aligned using ClustalW  69  , as implemented in the Geneious package  70. Internal Transcribed Spacer regions (ITS1 and ITS2) were identified and annotated based on a previous work on ITS secondary structures of MAST-4  53.
Gene prediction of co-assembled genomes and taxonomic profiling. The initial set of CEGs predicted with CEGMA were used to train the Augustus ab initio gene predictor  71  prior to its execution on the full co-assembly using defaults parameters (Supplementary Fig. S1d ). Genes were annotated using BLASTp (e-value < 10 -5 ) and rpsBLAST against, respectively, KEGG Orthology 54 and the Conserved Domains and Protein Classification database of NCBI (NCBI-CDD)  55. BLASTp hits with at least 100 bp alignments including at least 30% of query coverage and > 25% similarity were kept.
Acknowledgements This work was supported by the  US NSF  grants  DEB-1031049  and  OCE-821374  (to M. E. S. ), by the  ANR French projects Oceanomics  ( ANR-11-BTBR-0008 , to C. V. ),  France Génomique  ( ANR-10-INBS-09 , to P. W. ), and  Prometheus  ( ANR-09-PCS-GENM_217 , to O. J. ), by the  EU  project  SINGEK  ( H2020-MSCA-ITN-2015-675752 , to R. M. ), and by the  Spanish  project  MEFISTO  ( CTM2013-43767-P , MINECO). J. -F. M. was supported by a  Marie Curie Intra-European Fellowship  ( PIEF-GA-2012-331190 , EU). R. L. was supported by  Juan de la Cierva  ( JCI-2010-06594 ,  MINECO)  and  Ramón y Cajal fellowships  ( RYC-2013-12554 , MINECO). Computing resources were obtained through the  MARBITS  platform at the  ICM-CSIC  as well as through the  Red Española de Supercomputación. We appreciate the efforts of the  Single Cell Genomic Center of Bigelow  (https://scgc.