Microbiota processing in mothur, UCR Workshop 2018
Kim Dill-McFarland and Madison Cox
January 15 & 18, 2018
Updated January 18, 2018
Online version available at http://rpubs.com/maddieSC/mothur_SOP_UCR_Jan_2018
Madison’s email: mcox4@wisc.edu Carlos’ (CeNAT) email: cgamboa@cenat.ac.cr
Slides: https://drive.google.com/file/d/1z8AMMVraj4_Zb5dzcJUMQ5w_Ftdz7r7-/view?usp=sharing
Tips for this workshop
You can copy-paste from this document into the Terminal window. This will be faster and prevent typos. Depending on your computer and your ssh client, this may be done by
- Ctrl+v
- Apple+v
- Right-click on your mouse
- GREEN boxes contain code that you can copy and paste to run on your machine.
- WHITE boxes contain sample output of this code, and nothing will happen if you try to copy it into your console.
- ORANGE boxes contain commands that you should NOT execute, either because they are unnecessary for today, or because they take too long to do during this workshop and have been done for you.
Puede copiar y pegar desde este documento en la ventana Terminal. Esto será más rápido y evitar errores tipográficos. Dependiendo de su computadora y su cliente ssh, esto puede hacerse por
- Ctrl + v
- Apple + v
- Haz clic con el botón derecho del mouse
- Las cajas VERDES contienen código que puede copiar y pegar para ejecutar en su máquina.
- Las cajac BLANCAS contienen ejemplos de salida de este código, y nada sucederá si intentas copiarlo en tu consola.
- Las cajas ANARANJADAS contienen comandos que NO debe ejecutar, ya sea porque son innecesarios para hoy, o porque toman demasiado tiempo para hacer durante este taller y se han hecho por usted.
Common issues
Issues during this workshop
If at any time during this workshop you fall behind, all of the intermediate files are available in /home/capacitaciones/bioinformatica-UCR/. You will see that the files created by each mothur command are listed after that command under “Output File Names’. These are the files you need in order to skip that step.
For example, if you did not complete align.seqs and want to skip ahead, you would need to transfer the .align and .align.report files to your directory. To do this, you would use
mothur > system(cp /home/capacitaciones/bioinformatica-UCR/example.trim.contigs.good.unique.align ~/MiSeq_SOP)
mothur > system(cp /home/capacitaciones/bioinformatica-UCR/example.trim.contigs.good.unique.align.report ~/MiSeq_SOP)Some of the larger files are also zipped.
example.trim.contigs.qual.zipexample.trim.contigs.good.unique.align.zipexample.trim.contigs.good.unique.good.align.zipexample.final.dist.zip
Therefore after transfer, you would also need to unzip them.
mothur > system(cp /home/capacitaciones/bioinformatica-UCR/example.final.dist.zip ~/MiSeq_SOP)
mothur > system(unzip ~/MiSeq_SOP/example.final.dist.zip)Common mistakes
Typos The absolute most common mistake when working in mothur is typos. The mothur file names get very long and somewhat repetitive so please double-check! This can even occur when copy-pasting from this document if you accidentally copy a leading or lagging space around the command.
fasta and count_table don’t match Another common error seen in mothur occurs when your .fasta and .count_table files do not match. This looks like
"[ERROR]: yourSequence is in fileA but not in fileB, please correct."This often occurs because you forgot to include your count file in a command or accidentally included the wrong one. You can see this when mothur picks a file to use when you don’t include one.
...
Using /home/mcox/MiSeq_SOP/example.trim.contigs.fasta as input file for the fasta parameter.
...You should make sure the file mothur picked is the correct one. And if you did provide this file in the command, check for typos!
Another reason a mismatch can occur is a multi-processor function failing. One or more processors may fail but the overall function appears to complete as the other processors finish. Thus, a file can be incomplete which would cause a mismatch. Running out of RAM will cause this to happen. The only way to correct this is to re-run the failed process (usually align.seqs, pre.cluster, dist.seqs, or cluster.split) on fewer processors.
Problemas comunes Problemas durante este taller Si en cualquier momento durante este taller se queda atrás, todos los archivos intermedios están disponibles en / mnt / fdsafdsafdsa / MOTHUR. Verá que los archivos creados por cada comando mothur se enumeran después de ese comando en “Nombres de archivo de salida”. Estos son los archivos que necesita para omitir ese paso.
Por ejemplo, si no completó align.seqs y quiere omitir el avance, deberá transferir los archivos.align y .align.report a su directorio. Para hacer esto, usarías
Algunos de los archivos más grandes también están comprimidos.
example.trim.contigs.qual.zipexample.trim.contigs.good.unique.align.zipexample.trim.contigs.good.unique.good.align.zipexample.final.dist.zip
Por lo tanto, después de la transferencia, también necesitarás descomprimirlos.
Errores comunes ** Typos ** El error más común absoluto cuando se trabaja en mothur es errores tipográficos. Los nombres de los archivos de mothur son muy largos y algo repetitivos, así que ¡por favor revisen! Esto incluso puede ocurrir al copiar y pegar desde este documento si accidentalmente copia un espacio inicial o rezagado alrededor del comando.
** fasta y count_table no coinciden ** Otro error común visto en mothur ocurre cuando sus archivos .fasta y.count_table no coinciden. Esto se ve como
Esto ocurre a menudo porque olvidó incluir su archivo de conteo en un comando o accidentalmente incluyó el incorrecto. Puedes ver esto cuando mothur toma un archivo para usar cuando no incluyes uno.
Debes asegurarte de que el archivo elegido por mothur sea el correcto. ¡Y si proporcionó este archivo en el comando, compruebe si hay errores tipográficos!
Otro motivo por el que puede producirse una falta de coincidencia es una falla de la función de multiprocesador. Uno o más procesadores pueden fallar pero la función general parece completarse cuando los otros procesadores finalizan. Por lo tanto, un archivo puede estar incompleto, lo que causaría una falta de coincidencia. Quedando sin RAM hará que esto suceda. La única forma de corregir esto es volver a ejecutar el proceso fallido (normalmente align.seqs,pre.cluster, dist.seqs, ocluster.split) en pocos procesadores.
Introduction
Goal
The goal of this tutorial is to demonstrate standard operating procedures (SOP) for processing amplicon data that are generated using Illumina’s MiSeq platform with paired end reads. This SOP is a combination of the Schloss lab and Suen lab (UW-Madison) procedures for 2x250bp reads of the V4 region of the 16S rRNA gene. There are a number of steps, as noted in this SOP, which require modification if you are using a different amplicon and/or read length.
Sources and citation
This SOP is based on Pat Schloss’s MiSeq SOP with some notes and modifications by Dr. Kim Dill-McFarland and Madison Cox. The original SOP can be found at http://www.mothur.org/wiki/MiSeq_SOP.
Published work using this SOP should cite https://www.ncbi.nlm.nih.gov/pubmed/23793624
— Kozich JJ, Westcott SL, Baxter NT, Highlander SK, Schloss PD. (2013): Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform. Applied and Environmental Microbiology 79(17): 5112-20. doi: 10.1128/AEM.01043-13.
Programs
We will complete all sequence processing in the program mothur (https://www.mothur.org/wiki/Main_Page). Details on all commands used here as well as others can be found at http://www.mothur.org/wiki/Category:Commands.
We will use the command-line version of this software on a remote Linux computer. Windows PC users will need to download an ssh client for connecting to the remote computer. We recommend the free, personal edition of MobaXterm (http://mobaxterm.mobatek.net). This software is a combination of an enhanced terminal window and ssh client.
In addition to the ssh client, Windows PC users will also need to download a file unzipping utility which is able to hangle tar/gzipped archives. I strongly recommend 7zip, which is free here.
Introducción Gol El objetivo de este tutorial es demostrar ** procedimientos operativos estándar (SOP) ** para procesar datos de amplicones que se generan utilizando la plataforma MiSeq de Illumina con lecturas de pares emparejados. Este SOP es una combinación de los procedimientos Schloss lab y Suen lab (UW-Madison) para lecturas de 2x250 pb de la región V4 del gen 16S rRNA. Hay una serie de pasos, como se señala en este SOP, que requieren modificaciones si está utilizando un amplicón diferente y / o longitud de lectura.
Fuentes y citas Este SOP se basa en el SOP MiSeq de Pat Schloss con algunas notas y modificaciones del Dr. Kim Dill-McFarland y Madison Cox. El SOP original se puede encontrar en http://www.mothur.org/wiki/MiSeq_SOP.
El trabajo publicado que usa este SOP debe citar https://www.ncbi.nlm.nih.gov/pubmed/23793624
— Kozich JJ, Westcott SL, Baxter NT, Highlander SK, Schloss PD. (2013): Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform. Applied and Environmental Microbiology 79(17): 5112-20. doi: 10.1128/AEM.01043-13.
Data in this SOP
Data description
The Suen lab is interested in understanding how rumen (gut) microbes contribute to nutrition and milk production in dairy cows. In this example data set, we are looking at two samples: one solid- (S) and one liquid-phase (L) sample of digesta from the rumen (largest stomach compartment) of a lactating dairy cow. Obviously this is an artifically small data set which we are using to demonstrate the mothur functions without having to wait several minutes to several hours for commands to run.
We will analyze these data in this tutorial using an OTU (Operational Taxonomic Unit) approach. This approach is the most commonly used in microbiota research. Other approaches such as phylotype and phylogenetic can be found in the online Schloss SOP.
The data in this SOP were generated using 16S V4 primers from Kozich et al. 2013. These primers allow you to sequence up to 1536 samples in parallel using only 80 primers (32+48) and obtain sequence reads of the V4 region of the 16S rRNA gene. Please consult the supplementary methods of that manuscript for more information. You can also see the latest Schloss wet-lab SOP for generating libraries similar to these (https://github.com/SchlossLab/MiSeq_WetLab_SOP).
Datos en este SOP Descripción de los datos El laboratorio de Suen está interesado en comprender cómo los microbios del rumen (intestino) contribuyen a la nutrición y la producción de leche en las vacas lecheras. En este conjunto de datos de ejemplo, estamos viendo dos muestras: una muestra de sólidos (S) y una de fase líquida (L) de digesta del rumen (compartimiento estomacal más grande) de una vaca lechera lactante. Obviamente, este es un conjunto de datos artificialmente pequeño que estamos utilizando para demostrar las funciones mothur sin tener que esperar varios minutos o varias horas para ejecutar los comandos.
Analizaremos estos datos en este tutorial utilizando un enfoque OTU (unidad taxonómica operativa). Este enfoque es el más comúnmente utilizado en la investigación de la microbiota. Otros enfoques, como el filotipo y el filogenético, se pueden encontrar en línea [Schloss SOP] (http://www.mothur.org/wiki/MiSeq_SOP).
Files
From the MiSeq Sequences come off the MiSeq as pairs of fastq.gz files with each pair representing the two paired reads per sample.
fastq files contain both the sequence data (fasta) and the quality score data (qual). gz is a form of compression to make the files smaller. If you aren’t getting these files off the sequencer, then you likely have the software parameters set incorrectly. The default MiSeq file names are sample.name_S##_L001_R#_001.fastq where
sample.nameis the name you give the machine when running your library (i.e. pool of all samples together)S##is the sample number given by the machine to ensure all file names are uniqueL001is lane 1 (MiSeqs only have one lane as opposed to HiSeqs)R#is the forward (R1) or reverse (R2) read for each sample001is the replicate number, which will always be001as the MiSeq requires unique sample names (again, as opposed to the HiSeq)
For mothur For this tutorial, you will need our 4 example samples’ fastq.gz files as well as two databases for alignment and classification. All of these have been pre-downloaded onto Kabré. However, if you need the databases in the future, they can be found through the mothur website. The database files include
- mothur-formatted Silva sequences, reference alignment, and taxonomy
- mothur-formatted GreenGenes sequences, reference alignment, and taxonomy
From mothur When mothur modifies your data, it creates a new file for each step. For example, if you start out with Example.fasta and run screen.seqs to remove poor-quality sequences, mothur will create Example.good.fasta. Thus, you can see all the steps you’ve run just by looking at the file names.
Note: It is generally easiest to use the “current” option for many of the commands since the file names get very long. Because this tutorial is meant to show people how to use mothur at a very nuts and bolts level, we will only selectively use the current option to demonstrate how it works. Generally, we will use the full file names for this tutorial.
Archivos ** Del MiSeq ** Las secuencias salen de MiSeq como pares de archivos fastq.gz, y cada par representa las dos lecturas emparejadas por muestra.
Los archivos fastq contienen tanto los datos de secuencia (fasta) como los de puntuación de calidad (qual). gz es una forma de compresión para hacer que los archivos sean más pequeños. Si no obtiene estos archivos del secuenciador, es probable que tenga los parámetros del software configurados incorrectamente. Los nombres de archivo predeterminados de MiSeq son sample.name_S ## _ L001_R # _001.fastq donde
sample.namees el nombre que le da a la máquina cuando ejecuta su biblioteca (* i.e. * conjunto de todas las muestras juntas)S ##es el número de muestra proporcionado por la máquina para garantizar que todos los nombres de archivo sean únicosL001es el carril 1 (MiSeqs solo tiene un carril en lugar de HiSeqs)R #es la lectura directa (R1) o inversa (R2) para cada muestra001es el número de réplica, que siempre será001ya que MiSeq requiere nombres de muestra únicos (de nuevo, a diferencia del HiSeq)
** Para mothur ** Para este tutorial, necesitará nuestros 4 ejemplos de ejemplos de archivos ’fastq.gzy dos ** bases de datos ** para su alineación y clasificación. Todos estos han sido descargados previamente enKabré`. Sin embargo, si necesita las bases de datos en el futuro, se pueden encontrar a través del sitio web de mothur. Los archivos de la base de datos incluyen
- secuencias Silva formateadas en mothur, alineación de referencia y taxonomía
- Secuencias de GreenGenes con formato mothur, alineación de referencia y taxonomía
** De mothur ** Cuando mothur modifica sus datos, crea un nuevo archivo para cada paso. Por ejemplo, si comienza con Example.fasta y ejecuta screen.seqs para eliminar secuencias de baja calidad, mothur crearáExample.good.fasta. Por lo tanto, puedes ver todos los pasos que has ejecutado con solo mirar los nombres de los archivos.
** Nota : generalmente es más fácil usar la opción “actual” para muchos de los comandos, ya que los nombres de los archivos son muy largos. Debido a que este tutorial está destinado a mostrar a las personas cómo usar mothur en un nivel muy básico, solo utilizaremos selectivamente la opción actual para demostrar cómo funciona. En general, utilizaremos los nombres completos de archivo para este tutorial **.
Get Started
Download and run mothur on your own machine
Computational requirement
Make sure that your machine is powerful enough for the data set you wish to run. What we are doing today is a very, very small data set, which should not pose any issues on a resonably new standard laptop. However, if you were to run a full data set on a low-end personal computer, you would likely run out of memory and crash it. Not to mention that many steps in mothur processing are very time-intensive. If you can, find a way to access a large server to run full data sets. My laptop has 16Gb of RAM and 4 processors. With these specifications, running on 2 processors, a small data set of ~25 samples (V4) takes about a day. A large data set of 400+ samples may take weeks. A V3-4 data set of any reasonable size would likely crash my machine (see the “large distance matrix” section below for an explanation). So make sure that whatever machine you use is up to the task before you start.
Download the latest verison of mothur
The latest mothur release can be found here. Go to the github page and download the latest version. At time of writing, it was version 1.39.5. Make sure to keep track of what version of mothur you are running, as it is important to report in any manuscript using data processed in mothur. Download the zip file appropriate for your operating system and place it on your desktop. You may choose to move it later, but for today we will leave it there. Unzip the file to a folder on the desktop, and rename the folder “mothur-1.39.5”. Inside you will find the mothur executable (mothur.exe), as well as other executables which certain parts of the mothur workflow rely on. Congratulations, you have downloaded the program.
Connect to the Kabré server
We will remotely connect to the Kabré cluster for these exercises.
Connect from a Mac
- Open a
Terminal(if you are not sure where that is simply search for the word Terminal with “Spotlight” tool i.e. the top right magnifying glass icon) - Within the Terminal, type the following command adapted with your own username (shown here for
mcox):ssh -Y mcox@cluster.cenat.ac.cr - The
-Yoption allows passing of any graphics from the cluster. The XQartz program should be installed for this to work (See https://www.xquartz.org) Kabréis the name of the cluster.
Connect from a Windows PC
There are 2 ways to connect to a remote computer:
Standard connection
- Start MobaXterm
- In the main window, type the ssh command containing your username and the computer name. For example for
mcox, you would write:ssh -Y mcox@cluster.cenat.ac.cr - The
-Yoption allows passing of any graphics from the cluster. The XQartz program should be installed for this to work (See https://www.xquartz.org) Kabréis the name of the cluster.
GUI method connection
- Start MobaXterm
- In the QuickConnect text box (upper left), enter the name of the computer you want to connect to: e.g.
@cluster.cenat.ac.cr - Press the
Returnbutton - In the next window, click on the
SSHicon (upper left). Click OK - At the “Login:” prompt that appears within the main window, enter the username provided e.g.
mcox - Enter the password provided (Note: nothing gets printed)
- You may be asked if you want to save the password. Click NO
- You should now be logged in the remote computer on the main window.
Pasting in MobaXterm
By default, the pasting utility in MobaXterm is Shift+Insert. You can change this by going into the “Settings” pulldown menu and opening “Keyboard Shortcuts”. Find “Paste” and use the “Edit keyboard shortcut” pulldown to select “Ctrl + V” to make it normal.
Empezar Descargue y ejecute mothur en su propia máquina Requisito computacional Asegúrese de que su máquina sea lo suficientemente potente para el conjunto de datos que desea ejecutar. Lo que estamos haciendo hoy es un conjunto de datos muy, muy pequeño, que no debería plantear ningún problema en una computadora portátil estándar nueva y resonable. Sin embargo, si tuviera que ejecutar un conjunto completo de datos en una computadora personal de gama baja, es probable que se quede sin memoria y se bloquee. Sin mencionar que muchos pasos en el procesamiento de mothur requieren mucho tiempo. Si puede, encuentre una forma de acceder a un servidor grande para ejecutar conjuntos completos de datos. Mi laptop tiene 16Gb de RAM y 4 procesadores. Con estas especificaciones, ejecutándose en 2 procesadores, un pequeño conjunto de datos de ~ 25 muestras (V4) toma alrededor de un día. Un gran conjunto de datos de más de 400 muestras puede llevar semanas. Es probable que un conjunto de datos V3-4 de cualquier tamaño razonable bloquee mi máquina (consulte la sección “matriz de distancia grande” a continuación para obtener una explicación). Así que asegúrese de que cualquier máquina que use esté a la altura de la tarea antes de comenzar.
Descargar la última versión de mothur El último lanzamiento de mothur se puede encontrar [aquí.] (Https://www.mothur.org/wiki/Download_mothur#Our_latest_Release_of_Mothur) Vaya a la página de github y descargue la última versión. Al momento de escribir, era la versión 1.39.5. Asegúrese de realizar un seguimiento de la versión de mothur que está ejecutando, ya que es importante informar en cualquier manuscrito con los datos procesados en mothur. Descargue el archivo zip apropiado para su sistema operativo y colóquelo en su escritorio. Puede elegir moverlo más tarde, pero hoy lo dejaremos allí. Descomprime el archivo en una carpeta en el escritorio y cambia el nombre de la carpeta “mothur-1.39.5”. En su interior encontrará el ejecutable de mothur (mothur.exe), así como otros ejecutables de los que dependen ciertas partes del flujo de trabajo de mothur. Felicidades, has descargado el programa.
Conectarse al servidor de Kabré Nos conectaremos remotamente al clúster de Kabré para estos ejercicios.
Conectar desde una Mac 1. Abra un Terminal (si no está seguro de dónde está, simplemente busque la palabra Terminal con la herramienta" Spotlight “, es decir, el icono de la lupa superior derecha) 2. Dentro de la Terminal, escriba el siguiente comando adaptado con ** su propio usuario ** (que se muestra aquí para mcox): ssh -Y mcox@cluster.cenat.ac.cr 3. La opción -Y permite pasar cualquier gráfico del clúster. El programa XQartz debe instalarse para que esto funcione (Ver https://www.xquartz.org) 4. Kabré es el nombre del clúster.
Connect desde una PC con Windows Hay 2 formas de conectarse a una computadora remota:
** Conexión estándar **
- Inicie MobaXterm
- En la ventana principal, escriba el comando ssh que contiene ** su nombre de usuario ** y el nombre de la computadora. Por ejemplo para
mcox, usted escribiría:ssh -Y mcox@cluster.cenat.ac.cr - La opción
-Ypermite pasar cualquier gráfico del clúster. El programa XQartz debe instalarse para que esto funcione (Ver https://www.xquartz.org) Kabrées el nombre del clúster.
** Conexión del método GUI **
- Inicie MobaXterm
- En el cuadro de texto ** QuickConnect ** (arriba a la izquierda), ingrese el nombre de la computadora a la que desea conectarse: * por ejemplo, *
@cluster.cenat.ac.cr - Presione el botón
Return - En la siguiente ventana, haga clic en el ícono
SSH(arriba a la izquierda). Haga clic en Aceptar - En el indicador de “Inicio de sesión:” que aparece en la ventana principal, ingrese el nombre de usuario provisto * por ej. *
mcox - Ingrese la contraseña provista (* Nota *: nada se imprime)
- Es posible que se le pregunte si desea guardar la contraseña. Click NO
- Ahora debe iniciar sesión en la computadora remota en la ventana principal.
** Pegado en MobaXterm **
Por defecto, la herramienta de pegar en MobaXterm es Shift + Insertar. Puede cambiar esto yendo al menú desplegable “Configuración” y abriendo “Atajos de teclado”. Encuentra “Pegar” y utiliza el menú desplegable “Editar atajo de teclado” para seleccionar “Ctrl + V” para hacerlo normal.
Home directory
Once you have connected to the Kabré computer, you “land” within your “Home” directory. For mcox, this would be:
/home/mcoxIf at any point you get lost within the labyrinth of directories, simply type cd to come back.
If you would like to back up by a single directory, type cd ..
To know where your terminal is “looking” at any moment, you can invoke the print working directory command: pwd
$ prompt
The $ in the Terminal is simply a way for the computer to notify you that it is ready to receive typed commands. ** Note: The $ is not typed when issuing commands. It is already there.**
For example, mcox would see the following Terminal prompt:
[mcox@login-0 ~]$When we first log in, we put on a login node in the Kabré system. In order to run programs, we need to move on to a different node. We will be on Zarate today.
$ ssh zarate-0bOR
$ ssh zarate-0aNow your prompt should reflect that we are on a “zarate” node rather than a “login” node.
Project folder: Stay organized
It is a general good practice to create a project folder or working directory in order to separate files from various projects. We will use the name MiSeq_SOP. Therefore, our first task is to create this directory within the most logical place: our “Home” directory.
First we can verify that we are indeed within the “Home” directory with the print working directory command pwd.
$ pwd/home/mcoxThen, we create the directory with mkdir
Create a project folder
$ mkdir MiSeq_SOPWe can verify that this directory exists by checking the content of our “Home” directory with the list command ls
$ lsMiSeq_SOPNext, move into your MiSeq_SOP folder. This will be done with the change directory command cd into MiSeq_SOP:
Change directory
$ cd MiSeq_SOPWe can check that this worked with the command print working directory pwd. For mcox, this would look like:
$ pwd/home/mcox/MiSeq_SOP
Directorio de inicio Una vez que se haya conectado a la computadora Kabré," aterrizará “dentro de su directorio” Inicio “. Para mcox, esto sería:
Si en algún momento te pierdes dentro del laberinto de directorios, simplemente escribe cd para volver.
Si desea realizar una copia de seguridad en un único directorio, escriba cd ..
Para saber dónde está “buscando” su terminal en cualquier momento, puede invocar el comando del directorio de trabajo de impresión: pwd
$ prompt El $ en el Terminal es simplemente una forma de que la computadora le notifique que está listo para recibir comandos escritos. ** * Nota *: El $ no se escribe al emitir comandos. Ya está allí. **
Por ejemplo, mcox vería el siguiente promptTerminal:
Cuando iniciamos sesión por primera vez, instalamos un nodo de inicio de sesión en el sistema Kabré. Para ejecutar programas, necesitamos pasar a un nodo diferente. Estaremos en Zárate hoy.
Ahora su mensaje debe reflejar que estamos en un nodo “zarate” en lugar de un nodo “iniciar sesión”.
Carpeta del proyecto: Manténgase organizado Es una buena práctica general crear una carpeta de proyecto o un directorio de trabajo para separar archivos de varios proyectos. Utilizaremos el nombre MiSeq_SOP. Por lo tanto, nuestra primera tarea es crear este directorio en el lugar más lógico: nuestro directorio “Inicio”.
Primero podemos verificar que estamos efectivamente dentro del directorio “Inicio” con el comando * print working directory * pwd.
Entonces, creamos el directorio con mkdir
** Crea una carpeta de proyecto **
Podemos verificar que este directorio exista comprobando el contenido de nuestro directorio “Inicio” con el comando de lista ls
Luego, vaya a su carpeta MiSeq_SOP. Esto se hará con el comando * change directory * cd en MiSeq_SOP:
** Cambiar directorio **
Podemos verificar que esto funcionó con el comando * print working directory * pwd. Para mcox, esto se vería así:
Get data files
In this section, we will retrieve the data files for our 6 samples. Remember that there are R1 and R2 files so we will end up with 12 fastq.gz files.
Note: If you wanted to download files from the Internet, you could use the command web get wget followed by the URL, or web address of the file to download. For example: wget https://www.mothur.org/w/images/d/d6/MiSeqSOPData.zip
Sequence fastq
We will retrieve fastq files for our 2 samples (4 files). We will copy cp a zipped folder containing the files that has already been downloaded to the local computer located in the directory: /home/capacitaciones/bioinformatica-UCR/. Then, you will unzip these files and move mv them into your MiSeq_SOP folder.
Copy data files into home directory
$ cp /home/capacitaciones/bioinformatica-UCR/fastq.zip ~/MiSeq_SOPNote: ~ is a representation of the home directory, /home/mcox with your username
Nothing will be printed to the window so check that the file has been moved.
$ lsfastq.zipNow, unzip the folder to pull out all the fastq.gz files. Remember to change the “mcox” to your username
Unzip data files
$ unzip /home/mcox/Miseq_SOP/fastq.zipThis will place all of the files in a folder named the same thing as the zipped object (i.e. fastq).
lsfastq/
fastq.zipSo we need to move the fastq.gz files from fastq/ into MiSeq_SOP/. We can move all of the files ending in fastq.gz at once with the wildcard *
$ mv /home/mcox/MiSeq_SOP/fastq/*fastq.gz ~/MiSeq_SOPYou can then check the content of the directory with the list ls command.
$ lsfastq/
fastq.zip
L62.1_S149_L001_R1_001.fastq.gz
L62.1_S149_L001_R2_001.fastq.gz
S62.1_S27_L001_R1_001.fastq.gz
S62.1_S27_L001_R2_001.fastq.gzFor these samples, the sample names are coded as sampleCow.day. For example, L62.1 is the liquid sample from cow #62 on day 1. These sequences are 250 bp and overlap in the V4 region of the 16S rRNA gene; this region is about 253 bp long.
Obtener archivos de datos En esta sección, recuperaremos los archivos de datos para nuestras 6 muestras. Recuerde que hay archivos R1 yR2, así que terminaremos con 12 archivos fastq.gz.
** Nota **: si desea descargar archivos de Internet, puede usar el comando * web get * wget seguido de la URL o la dirección web del archivo para descargar. Por ejemplo: wget https: // www.mothur.org / w / images / d / d6 / MiSeqSOPData.zip
Secuencia fastq Recuperaremos archivos fastq para nuestras 2 muestras (4 archivos). Copiaremos cp una carpeta comprimida que contiene los archivos que ya se han descargado a la computadora local ubicada en el directorio:/ mnt / software / workshop / data / MOTHUR /. Luego, descomprimirás estos archivos y moverás mv en tu carpetaMiSeq_SOP.
** Nota **: ~ es una representación del directorio de inicio, / home / mcox con su número de estudiante
No se imprimirá nada en la ventana, así que verifique que el archivo se haya movido.
Ahora, descomprime la carpeta para extraer todos los archivos fastq.gz. Recuerde ** cambiar el “mcox” a su nombre de usuario **
Esto colocará todos los archivos en una carpeta llamada lo mismo que el objeto comprimido (* i.e. * fastq).
Así que tenemos que mover los archivos fastq.gz defastq /aMiSeq_SOP /. Podemos mover todos los archivos que terminan en fastq.gz a la vez con el comodín*
Luego puede verificar el contenido del directorio con el comando list ls.
Para estas muestras, los nombres de las muestras se codifican como sampleCow.day. Por ejemplo, L62.1 es la muestra líquida de la vaca # 62 en el día 1. Estas secuencias son de 250 pb y se superponen en la región V4 del gen 16S rRNA; esta región tiene aproximadamente 253 pb de largo.
Databases
We will now copy the Silva and GreenGenes database files into your home folders with cp. For alignment, we will be using a reduced version of Silva that only contains the V4 region of 16S since we know that is our amplicon region. A full Silva database is available at https://www.mothur.org/wiki/Silva_reference_files. Also, because silva.v4.align is large, it is zipped and will need to be unzipped once in your folder.
All of these files are located in /opt/src/SilvaDB and /opt/src/GreegenesDB
- silva.v4.align.zip
- silva.nr_v128.tax
- gg_13_8_99.fasta
- gg_13_8_99.gg.tax
Retrieve databases
$ cp /opt/src/SilvaDB/silva.v128.v4.align.zip ~/MiSeq_SOP
$ cp /opt/src/SilvaDB/silva.nr_v128.tax ~/MiSeq_SOP
$ cp /opt/src/GreegenesDB/gg_13_8_99.fasta ~/MiSeq_SOP
$ cp /opt/src/GreegenesDB/gg_13_8_99.gg.tax ~/MiSeq_SOPUnzip align file
$ unzip ~/MiSeq_SOP/silva.v128.v4.align.zipsilva.v128.v4.align.zip was a single zipped file and therefore, was not placed in a new folder like the fastq.gz files. So, it does not need to be moved.
If you do a final check of your folder, you should have all of the following files.
$ lsfastq/
fastq.zip
gg_13_8_99.fasta
gg_13_8_99.gg.tax
L62.1_S149_L001_R1_001.fastq.gz
L62.1_S149_L001_R2_001.fastq.gz
S62.1_S27_L001_R1_001.fastq.gz
S62.1_S27_L001_R2_001.fastq.gz
silva.nr_v128.tax
silva.v4.align
silva.v4.align.zip
Bases de datos Ahora copiaremos los archivos de la base de datos de Silva y GreenGenes en sus carpetas de inicio con cp. Para la alineación, utilizaremos una versión reducida de Silva que solo contiene la región V4 de 16S, ya que sabemos que es nuestra región de amplicón. Una base de datos completa de Silva está disponible en https://www.mothur.org/wiki/Silva_reference_files. Además, como silva.v4.align es grande, está comprimido y deberá descomprimirse una vez en su carpeta.
Todos estos archivos se encuentran en /opt/src/SilvaDB y /opt/src/GreegenesDB
- silva.v4.align.zip
- silva.nr_v128.tax
- gg_13_8_99.fasta
- gg_13_8_99.gg.tax
** Recuperar bases de datos **
** Descomprima el archivo de alineación **
silva.v128.v4.align.zip era un único archivo comprimido y, por lo tanto, no se colocó en una carpeta nueva como los archivosfastq.gz. Por lo tanto, no necesita ser movido.
Si realiza una verificación final de su carpeta, debe tener todos los siguientes archivos.
Begin in mothur
The software mothur will be accessed via command-line on the Terminal. We can check that the computer knows about the software with the command which. This will echo back the location of the software:
$ module load mothur/1.39.5
$ which mothur/opt/bioinf/mothur-1.39.5/mothurWe can also verify the version that is installed. It is important to note the version used as the calculations or command syntax may differ slightly with different versions.
We can check what version is installed with the following command that will print out the results on the screen:
$ mothur --versionLinux 64Bit Version
Mothur version=1.39.5
Release Date=3/20/2017
Empieza en mothur El software mothur se accederá a través de la línea de comandos en la Terminal. Podemos verificar que la computadora conozca el software con el comando which. Esto recordará la ubicación del software:
También podemos verificar la versión que está instalada. Es importante tener en cuenta la versión utilizada ya que los cálculos o la sintaxis del comando pueden diferir ligeramente con las diferentes versiones.
Podemos verificar qué versión está instalada con el siguiente comando que imprimirá los resultados en la pantalla:
General mothur syntax
mothur automatically creates file names based on commands run with that file.
Sometimes new versions slightly change the output names which cause later commands to fail as they cannot find the correct file. For example, version 1.35 outputs the file ...precluster.uchime.accnos after chimera.uchime while version 1.36 outputs ...precluster.denovo.uchime.accnos because new options were added and it now specifies that you used the denovo method.
If you are having issues, the first thing to check is whether or not all the needed input files exist and are named exactly what they are in your command. If they are not, make the necessary changes to the next command and try re-running.
We will run these exercises on the Kabré cluster. These resources have to be shared among all the attendees in class. I will let you know what number of processors we will be using for today before we begin.
If you were to run on your desktop or laptop computer on your own the following recommendations from Kim would be useful:
— Before you begin, there are some things to consider. If this is your first time using mothur, I recommend running through this tutorial will only 2 samples so that you can get used to command line and identify any problem areas without having to wait for computationally intensive steps to run. This will also minimize RAM load so that you can practice on your own computer before moving to larger server usage (more about that below).
— Once you move up to a bigger data set, you need to consider the computer you’re running mothur on. Many commands can be run on multiple processors by adding processors = #. However, you can only use as many as your computer actually has. A good rule of thumb is to not run mothur using more than (all of your processors - 2). This leaves 1 processor for your OS and 1 for any other programs you may have open. Try to minimize other things open if possible. The other thing to consider is RAM. If you have very little RAM, only run mothur with 1 processor to minimize RAM load. This will make this very slow going and a laptop might still run out of RAM and crash. It’s impossible to say exactly how much RAM you will need for a given data set. So, when thinking about running a large data set consisting of an entire MiSeq run or more, consider getting access to 0.5TB+ RAM machines. Madison has these resources available through the Center for High Throughput Computing, http://chtc.cs.wisc.edu/.
Sintaxis general de mothur mothur crea automáticamente nombres de archivo basados en comandos ejecutados con ese archivo.
A veces, las versiones nuevas cambian ligeramente los nombres de salida que hacen que los comandos posteriores fallen ya que no pueden encontrar el archivo correcto. Por ejemplo, la versión 1.35 genera el archivo ... precluster.uchime.accnos después de chimera.uchime, mientras que la versión 1.36 muestra... precluster.denovo.uchime.accnos porque se agregaron nuevas opciones y ahora especifica que usted utilizó el método denovo
Si tiene problemas, lo primero que debe comprobar es si existen o no todos los archivos de entrada necesarios y cómo se llaman exactamente en el comando. Si no lo son, realice los cambios necesarios en el próximo comando e intente volver a ejecutarlo.
Realizaremos estos ejercicios en el clúster Kabré. Estos recursos deben ser compartidos entre todos los asistentes en clase. Le informaré la cantidad de procesadores que utilizaremos hoy antes de comenzar.
Si tuviera que ejecutar en su computadora de escritorio o portátil por su cuenta, las siguientes recomendaciones de Kim serían útiles:
— Antes de comenzar, hay algunas cosas que considerar. Si esta es la primera vez que usa mothur, le recomiendo que ejecute este tutorial solo con 2 muestras para que pueda acostumbrarse a la línea de comandos e identificar cualquier área problemática sin tener que esperar a que se ejecuten los pasos computacionalmente intensivos. Esto también reducirá al mínimo la carga de RAM para que pueda practicar en su propia computadora antes de pasar a un uso más grande del servidor (más sobre eso más adelante).
— Una vez que pasas a un conjunto de datos más grande, debes considerar la computadora con la que estás ejecutando Mothur. Muchos comandos se pueden ejecutar en procesadores múltiples al agregar procesadores = #. Sin embargo, solo puedes usar tantos como tu computadora realmente tenga. Una buena regla general es no ejecutar mothur usando más de (todos sus procesadores - 2). Esto deja 1 procesador para su sistema operativo y 1 para cualquier otro programa que pueda tener abierto. Intenta minimizar otras cosas abiertas si es posible. La otra cosa a considerar es RAM. Si tiene muy poca RAM, solo ejecute mothur con 1 procesador para minimizar la carga de RAM. Esto hará que esto sea muy lento y una computadora portátil aún se quede sin memoria RAM y se bloquee. Es imposible decir exactamente la cantidad de RAM que necesitará para un conjunto de datos determinado. Por lo tanto, cuando piense en ejecutar un conjunto de datos grande que consista en una ejecución completa de MiSeq o más, considere obtener acceso a máquinas de 0.5TB + RAM. Madison tiene estos recursos disponibles a través del Centro de Computación de Alto Rendimiento, http://chtc.cs.wisc.edu/.
Log files
Every data set is different, and every data set presents its own set of problems. Logfiles are automatically generated by mothur and help to determine where issues may have arisen if running in batch mode. They also help with recordkeeping and reproducibility. Your log file is a live document which is continually updated over the course of a session. If you close mothur and reopen, an new logfile will be created. The logs are only differentiated by a long string of numbers, like mothur.1490885401.logfile so it may be a good idea to rename logfiles to be more informative, such as indicating dates and data sets.
Archivos de registro Cada conjunto de datos es diferente, y cada conjunto de datos presenta su propio conjunto de problemas. Loghuriles son generados automáticamente por mothur y ayudan a determinar dónde pueden haber surgido problemas si se ejecuta en modo batch. También ayudan con el mantenimiento de registros y la reproducibilidad. Su archivo de registro es un documento en vivo que se actualiza continuamente durante el transcurso de una sesión. Si cierra mothur y vuelve a abrir, se creará un nuevo archivo de registro. Los registros solo se diferencian por una larga cadena de números, como mothur.1490885401.logfile, por lo que puede ser una buena idea cambiar el nombre de los archivos de registro para que sean más informativos, como indicar fechas y conjuntos de datos.
Start mothur via command line
We can now start the program. We will all call the same path to execute mothur in our own folders. Since we are each in our own MiSeq_SOP folders, mothur will automatically save your files to your folder. Note that once this command is issued, we will be working within the mothur software as reminded by the new prompt mothur > which will only accept mothur line commands.
Start mothur
$ /opt/bioinf/mothur-1.39.5/mothurmothur v.1.39.5
Last updated: 3/20/2017
by
Patrick D. Schloss
Department of Microbiology & Immunology
University of Michigan
http://www.mothur.org
When using, please cite:
Schloss, P.D., et al., Introducing mothur: Open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol, 2009. 75(23):7537-41.
Distributed under the GNU General Public License
Type 'help()' for information on the commands that are available
For questions and analysis support, please visit our forum at https://www.mothur.org/forum
Type 'quit()' to exit program
mothur >
Inicie mothur a través de la línea de comando Ahora podemos comenzar el programa. Todos llamaremos el mismo camino para ejecutar mothur en nuestras propias carpetas. Como cada uno de nosotros está en nuestras propias carpetas MiSeq_SOP, mothur guardará automáticamente sus archivos en su carpeta. Tenga en cuenta que una vez que se emite este comando, trabajaremos dentro del software mothur tal como lo recuerda el nuevo mensaje mothur>, que solo aceptará comandos de línea mothur.
** Iniciar mothur **
Basic sequence clean-up
Combine paired end reads
Paired-end runs on the Illumina MiSeq yield two reads for every sequence. One running in from the 5 prime end (read 1, R1) and one from the 3 prime end (read 2, R2). You need to combine these reads so that there is a single, longer read for every sequence (a contig). You do this by determining the overlapping regions and combining with the make.contigs command.
This command will: * extract the sequence and quality score data from your ‘fastq’ files * create the reverse complement of the reverse read and then * join the reads into contigs
We have a very simple algorithm to do this. First, we align the pairs of sequences. Next, we look across the alignment and identify any positions where the two reads disagree. If one sequence has a base and the other has a gap, the quality score of the base must be over 25 to be considered real. If both sequences have a base at that position, then we require one of the bases to have a quality score 6 or more points better than the other. If it is less than 6 points better, then we set the consensus base to an N.
You could run make.contigs on each sample individually and then combine the outputs. However, this is very time consuming. So, make.contigs can also take a tab-delimited file listing all the samples and their corresponding file names. We will have mothur make this file using make.file. Remember to change to your username!
Importantly, the make.contigs command CAN accept the gzipped files which are generated by the MiSeq IF your computer is able to unzip them. So, if you are running mothur on OSX or Linux, you will have no issues as long as you let mothur know to look for the gzipped filetype. However, if you are running mothur on Windows, you will encounter all sorts of issues using gzipped files. You will have to use some sort of utility, like 7zip (http://www.7-zip.org/), to extract the fastq files prior to beginning. In this case, you would omit the “type=gz” in the following command, and mothur will look for the default “.fastq” filetype. Because Kabré is a Linux machine, we will be able to use the gzipped files with no issue in thei workshop.
make.file
mothur > make.file(inputdir=/home/mcox/MiSeq_SOP, type=gz)Setting input directory to: /home/mcox/MiSeq_SOP/
Output File Names:
/home/mcox/MiSeq_SOP/stability.filesThis will output a .files which is setup as
sample name -tab- R1file -tab- R2file
Let’s take a look at it.
mothur > system(head -5 stability.files)MiSeq /home/mcox/MiSeq_SOP/L62.1_S149_L001_R1_001.fastq.gz /home/mcox/MiSeq_SOP/L62.1_S149_L001_R2_001.fastq.gz
MiSeq /home/mcox/MiSeq_SOP/S62.1_S27_L001_R1_001.fastq.gz /home/mcox/MiSeq_SOP/S62.1_S27_L001_R2_001.fastq.gzNote: Notice that if you want to perform a terminal command outside of mothur while the program is running, you can use the system() command
Note: If you name two sets of files the same, they will be merged into one sample during make.contigs. This is useful if you sequenced a sample more than once.
We can see here that as a default in this version of mothur, both of our samples have been named “MiSeq”. This has to do with a change in the labeling of the .fastq files which has caused mothur to misinterperet sample names. In a normal case, you would edit the file in a text editor to give the correct sample names in the first column, and prevent ALL of your samples from grouping together under a single sample named “MiSeq”. For today, I have created a corrected file for you to use so that we can move forward. We are going to copy this file over from the same place we grabbed the fastq files.
Copy fixed .files file, look at it again
mothur > system(cp /home/capacitaciones/bioinformatica-UCR/stability_fixed.files ~/MiSeq_SOP)
mothur > system(head -5 stability_fixed.files)L62 L62.1_S149_L001_R1_001.fastq.gz L62.1_S149_L001_R2_001.fastq.gz
S62 S62.1_S27_L001_R1_001.fastq.gz S62.1_S27_L001_R2_001.fastq.gzNow we have a correctly formatted file, with the first column containing unique sample names, the second containing the filename (or full file path) for our forward reads, and the third column containing the filename (or full file path) for our reverse reads.
Note: Your sample names cannot contain a dash (-) as mothur uses this symbol to denote lists in some commands.
The default file name of “stability” is not a useful name for us so we’re going to re-name this file so that all subsequent files are named correctly.
If you are working on your own machine, rather than on a server, you can make this change with a regular text editor. I like to use Excel.
Rename ‘.files’
mothur > system(cp stability_fixed.files example.files)Note: If you were running on a Windows PC, you would use copy instead of cp since the system command calls upon the operating system of the computer. For this SOP, we are all running on a Linux server
Now, let’s begin our data clean-up by combining R1 and R2 with make.contigs. Please wait for instructions on how many processors to use as this will vary depending on the number of people attending today.
Limpieza de secuencia básica Combine las lecturas finales pareadas Las ejecuciones de extremos emparejados en Illumina MiSeq producen dos lecturas para cada secuencia. Uno que se ejecuta desde el 5 principal (leer 1, ‘R1’) y uno desde el 3 principal (leer 2, ‘R2’). Debe combinar estas lecturas para que haya una sola lectura más larga para cada secuencia (un contig). Para ello, determine las regiones superpuestas y combine con el comando make.contigs.
! [] (pairs_01.png)
Este comando: * extrae la secuencia y los datos de puntaje de calidad de tus archivos ‘fastq’ * crea el complemento inverso de la lectura inversa y luego * une las lecturas en contigs
Tenemos un algoritmo muy simple para hacer esto. Primero, alineamos los pares de secuencias. Luego, observamos la alineación e identificamos cualquier posición en la que las dos lecturas no concuerden. Si una secuencia tiene una base y la otra tiene una brecha, la puntuación de calidad de la base debe ser superior a 25 para que se considere real. Si ambas secuencias tienen una base en esa posición, entonces requerimos que una de las bases tenga un puntaje de calidad de 6 o más puntos mejor que el otro. Si es menos de 6 puntos mejor, establecemos la base de consenso en N.
Puede ejecutar make.contigs en cada muestra individualmente y luego combinar las salidas. Sin embargo, esto consume mucho tiempo. Por lo tanto, make.contigs también puede tomar un archivo delimitado por tabulaciones que enumera todas las muestras y sus nombres de archivo correspondientes. Tendremos que mothur haga este archivo usando make.file. ** ¡Recuerde cambiar su número de estudiante! **
** Es importante destacar que ** el comando make.contigs PUEDE aceptar los archivos gzip que genera MiSeq SI su computadora puede descomprimirlos. Por lo tanto, si está ejecutando mothur en OSX o Linux, no tendrá problemas siempre y cuando deje que mothur sepa que debe buscar el tipo de archivo gzip. Sin embargo, si está ejecutando mothur en Windows, encontrará todo tipo de problemas al usar archivos gzip. Deberá utilizar algún tipo de utilidad, como 7zip (http://www.7-zip.org/), para extraer los archivos fastq antes de comenzar. En este caso, omitiría el “tipo = gz” en el siguiente comando, y mothur buscará el tipo de archivo “.fastq” predeterminado. Debido a que Kabré es una máquina Linux, podremos usar los archivos comprimidos sin problemas en su taller.
** make.file **
Esto generará un .files que se configura como
nombre de la muestra -tab- R1file -tab-R2file
Echemos un vistazo.
** Nota : tenga en cuenta que si desea ejecutar un comando de terminal fuera de mothur mientras el programa se está ejecutando, puede usar el comando system () **
** Nota **: si nombra dos conjuntos de archivos iguales, se fusionarán en una muestra durante make.contigs. Esto es útil si secuencia una muestra más de una vez.
Aquí podemos ver que, como valor predeterminado en esta versión de mothur, nuestras dos muestras se han denominado “MiSeq”. Esto tiene que ver con un cambio en el etiquetado de los archivos .fastq que ha provocado que mothur malinterprete nombres de muestra. En un caso normal, debe editar el archivo en un editor de texto para dar los nombres de muestra correctos en la primera columna y evitar que TODAS sus muestras se agrupen en una sola muestra llamada “MiSeq”. Para hoy, he creado un archivo corregido para que lo use y podamos avanzar. Vamos a copiar este archivo desde el mismo lugar donde tomamos los archivos fastq.
** Copia el archivo .files fijo, míralo de nuevo **
Ahora tenemos un archivo correctamente formateado, con la primera columna que contiene nombres de muestra únicos, el segundo contiene el nombre de archivo (o la ruta completa) para nuestras lecturas hacia adelante y la tercera columna contiene el nombre de archivo (o la ruta completa) para nuestras lecturas inversas .
** Nota **: Sus nombres de muestra no pueden contener un guion (-) ya que mothur usa este símbolo para indicar listas en algunos comandos.
El nombre de archivo predeterminado de “estabilidad” no es un nombre útil para nosotros, por lo que vamos a renombrar este archivo para que todos los archivos posteriores tengan el nombre correcto.
Si está trabajando en su propia máquina, en lugar de en un servidor, puede hacer este cambio con un editor de texto normal. Me gusta usar Excel.
** Renombrar ‘archivos’ **
** Nota **: si estuviera ejecutando en una PC con Windows, usaría copy en lugar decp ya que el comando del sistema invoca el sistema operativo de la computadora. Para este SOP, todos estamos corriendo en un servidor Linux
Ahora, comencemos nuestra limpieza de datos combinando R1 yR2 con make.contigs. ** Espere las instrucciones ** sobre cuántos procesadores usar, ya que esto variará según la cantidad de personas que asisten hoy.
make.contigs
mothur > make.contigs(file=example.files, processors=2)
Using 2 processors.
>>>>> Processing file pair /home/mcox/MiSeq_SOP/L62.1_S149_L001_R1_001.fastq.gz - /home/mcox/MiSeq_SOP/L62.1_S149_L001_R2_001.fastq.gz (files 1 of 2) <<<<<
>>>>> Processing file pair /home/mcox/MiSeq_SOP/S62.1_S27_L001_R1_001.fastq.gz - /home/mcox/MiSeq_SOP/S62.1_S27_L001_R2_001.fastq.gz (files 2 of 2) <<<<<
Making contigs...
Making contigs...
1000
[truncated]
It took 103 secs to assemble 130107 reads.
It took 107 secs to process 145438 sequences.
Group count:
L62 15331
S62 130107
Total of all groups is 145438
Output File Names:
/home/mcox/MiSeq_SOP/example.trim.contigs.fasta
/home/mcox/MiSeq_SOP/example.trim.contigs.qual
/home/mcox/MiSeq_SOP/example.contigs.report
/home/mcox/MiSeq_SOP/example.scrap.contigs.fasta
/home/mcox/MiSeq_SOP/example.scrap.contigs.qual
/home/mcox/MiSeq_SOP/example.contigs.groups
[WARNING]: your sequence names contained ':'. I changed them to '_' to avoid problems in your downstream analysis.The Warning simply explains that the sequence names have been altered by changing the colon character into the underscore character because the colon character would cause problems for phylogenetic software, as would the dash character - already explained above.
This command will also produce several files that you will need down the road: example.trim.contigs.fasta and example.contigs.groups. These contain the sequence data and group identity for each sequence.
The example.contigs.report file will tell you something about the contig assembly for each read or why it failed. The fileexample.scrap.contigs.fasta contains sequences that fail to make contigs. The file example.contigs.groups contains a list of sequence names (given by the machine) and which sample they belong to.
We should now take a look at our data. This is done repeatedly in this SOP using ‘summary.seqs’
Note: mothur is smart enough to remember that we used 3 processors in make.contigs and so it will use that throughout your current session unless you specify otherwise.
summary.seqs
mothur > summary.seqs(fasta=example.trim.contigs.fasta)Using 2 processors.
Start End NBases Ambigs Polymer NumSeqs
Minimum: 1 248 248 0 3 1
2.5%-tile: 1 252 252 0 3 3636
25%-tile: 1 253 253 0 3 36360
Median: 1 253 253 0 4 72720
75%-tile: 1 253 253 0 4 109079
97.5%-tile: 1 254 254 5 5 141803
Maximum: 1 495 495 81 126 145438
Mean: 1 252.935 252.935 0.441852 3.83601
# of Seqs: 145438
Output File Names:
/home/mcox/MiSeq_SOP/example.trim.contigs.summary
It took 1 secs to summarize 145438 sequences.This output tells us that we have 145,438 sequences that are ~253 bases long.
Things to look out for
- The longest read in the data set is 495 bp. Recall that the reads are supposed to be ~250 bp each. This read clearly didn’t assemble well (or at all).
- ~2.5% of our sequences had some ambiguous base calls (i.e. N).
- ~2.5% of our sequences have long homopolymers (i.e. AAAAAAAAA)
We’ll take care of these issues in the next step when we run screen.seqs.
** make.contigs **
La * Advertencia * simplemente explica que los nombres de las secuencias han sido alterados cambiando el carácter de dos puntos al carácter de guión bajo porque el carácter de dos puntos causaría problemas para el software filogenético, como lo haría el personaje del tablero, ya explicado anteriormente.
Este comando también producirá varios archivos que necesitará en el futuro: example.trim.contigs.fasta yexample.contigs.groups. Estos contienen los datos de secuencia y la identidad de grupo para cada secuencia.
El archivo example.contigs.report le dirá algo sobre el conjunto de contig para cada lectura o por qué falló. El archivo example.scrap.contigs.fasta contiene secuencias que no pueden hacer contigs. El archivo example.contigs.groups contiene una lista de nombres de secuencia (dados por la máquina) y a qué muestra pertenecen.
Ahora deberíamos echar un vistazo a nuestros datos. Esto se hace de forma repetida en este SOP usando ‘summary.seqs’
** Nota **: mothur es lo suficientemente inteligente como para recordar que utilizamos 3 procesadores en make.contigs y así lo usará durante toda la sesión actual a menos que especifique lo contrario.
** summary.seqs **
Este resultado nos dice que tenemos 145,438 secuencias que tienen ~ 253 bases de longitud.
Cosas a tener en cuenta
- La lectura más larga en el conjunto de datos es de 495 pb. Recuerde que las lecturas son ~ 250 bp cada una. Esta lectura claramente no se ensambló bien (o en absoluto).
- ~ 2.5% de nuestras secuencias tenían algunas llamadas base ambiguas (* i.e. * N).
- ~ 2.5% de nuestras secuencias tienen homopolímeros largos (* i.e. * AAAAAAAAA)
Nos ocuparemos de estos problemas en el siguiente paso cuando ejecutamos screen.seqs.
Remove poor quality
This implementation of the command will remove sequences with any ambiguous bases, anything longer than 300 bp, and anything with a homopolymer longer than 8 bases.
screen.seqs
mothur > screen.seqs(fasta=example.trim.contigs.fasta, group=example.contigs.groups, maxambig=0, maxlength=300, maxhomop=8)Processing sequence: 100
Processing sequence: 200
[truncated]
Processing sequence: 72700
Processing sequence: 72773
Output File Names:
/home/mcox/MiSeq_SOP/example.trim.contigs.good.fasta
/home/mcox/MiSeq_SOP/example.trim.contigs.bad.accnos
/home/mcox/MiSeq_SOP/example.contigs.good.groups
It took 3 secs to screen 145438 sequences.Alternate command: There’s another way to run this using the output from summary.seqs. This may be faster because the summary.seqs output file already has the number of ambiguous bases and sequence length calculated for your sequences.
mothur > screen.seqs(fasta=example.trim.contigs.fasta, group=example.contigs.groups, summary=example.trim.contigs.summary, maxambig=0, maxlength=300, maxhomop=8)
Eliminar la mala calidad Esta implementación del comando eliminará las secuencias con bases ambiguas, cualquier cosa más larga que 300 pb, y cualquier cosa con un homopolímero de más de 8 bases.
** screen.seqs **
Comando alternativo: hay otra forma de ejecutar esto usando la salida de summary.seqs. Esto puede ser más rápido porque el archivo de salida summary.seqs ya tiene el número de bases ambiguas y la longitud de secuencia calculada para sus secuencias.
Using “current”
Now that we’ve created some files in mothur, we can ask mothur what it knows about our current session.
get.current
mothur > get.current()Current RAM usage: 0.0287514 Gigabytes. Total Ram: 15.5728 Gigabytes.
Current files saved by mothur:
fasta=/home/mcox/MiSeq_SOP/example.trim.contigs.good.fasta
group=/home/mcox/MiSeq_SOP/example.contigs.good.groups
qfile=/home/mcox/MiSeq_SOP/example.trim.contigs.qual
contigsreport=/home/mcox/MiSeq_SOP/example.contigs.report
processors=2
summary=/home/mcox/MiSeq_SOP/example.trim.contigs.summary
file=/home/mcox/MiSeq_SOP/stability.files
Current input directory saved by mothur: /home/mcox/MiSeq_SOP/
Current default directory saved by mothur: /opt/bioinf/mothur-1.39.5/
Current working directory: /home/mcox/MiSeq_SOP/
Output File Names:
current_files.summaryWhat this means is that mothur remembers your latest fasta and group file as well as the number of processors you are running on. This is useful because you can shorten commands by telling mothur to use the current files. So all of the following are equivalent at this point in analysis and you would get the same output for each.
mothur > summary.seqs(fasta=example.trim.contigs.good.fasta)
mothur > summary.seqs(fasta=current)
mothur > summary.seqs()Using 2 processors.
Start End NBases Ambigs Polymer NumSeqs
Minimum: 1 250 250 0 3 1
2.5%-tile: 1 252 252 0 3 3139
25%-tile: 1 253 253 0 3 31386
Median: 1 253 253 0 4 62771
75%-tile: 1 253 253 0 4 94156
97.5%-tile: 1 254 254 0 5 122402
Maximum: 1 294 294 0 8 125540
Mean: 1 252.889 252.889 0 3.82739
# of Seqs: 125540
Output File Names:
/home/mcox/MiSeq_SOP/example.trim.contigs.good.summary
It took 2 secs to summarize 125540 sequences.Note: If you exit out of mothur and then re-start, it will not remember the current files from your last session.
For the purposes of this tutorial, we will write out the names of the files for clarity.
Uso de “current” Ahora que hemos creado algunos archivos en mothur, podemos preguntarle a Mothur qué sabe sobre nuestra sesión actual.
** get.current **
Lo que esto significa es que mothur recuerda su último archivo fasta ygroup, así como la cantidad de procesadores en los que se está ejecutando. Esto es útil porque puede acortar los comandos al decirle a mothur que use los archivos actuales. Entonces, todo lo siguiente es equivalente en este punto del análisis y obtendría el mismo resultado para cada uno.
** Nota **: si sale de mothur y luego vuelve a comenzar, no recordará los archivos actuales de su última sesión.
A los fines de este tutorial, escribiremos los nombres de los archivos para mayor claridad.
Unique sequences
At this point, our sequencing error rate has probably dropped more than an order of magnitude and we have 125,540 sequences We anticipate that many of our sequences are duplicates of each other. Because it’s computationally wasteful to process the same thing a bazillion times, we’ll select only unique sequences using the unique.seqs command. This way, the computer can analyze one sequence and then apply that outcome to every identical sequence in the data set.
unique.seqs
mothur > unique.seqs(fasta=example.trim.contigs.good.fasta)[truncated]
125000 20396
125540 20469
Output File Names:
/home/mcox/MiSeq_SOP/example.trim.contigs.good.names
/home/mcox/MiSeq_SOP/example.trim.contigs.good.unique.fastaIf two sequences have the same identical sequence, then they’re considered duplicates and will get merged.
In the screen output there are two columns:
- the first is the number of sequences characterized and
- the second is the number of unique sequences remaining.
So after running unique.seqs we have gone from 125,540 total to 20,469 unique sequences. This will make our life much easier moving forward.
Secuencias únicas En este punto, nuestra tasa de error de secuenciación probablemente ha disminuido más de un orden de magnitud y tenemos 125,540 secuencias. Anticipamos que muchas de nuestras secuencias son duplicadas entre sí. Debido a que es un desperdicio computacional procesar lo mismo miles de veces, seleccionaremos solo secuencias únicas usando el comando unique.seqs. De esta forma, la computadora puede analizar una secuencia y luego aplicar ese resultado a cada secuencia idéntica en el conjunto de datos.
** unique.seqs **
Si dos secuencias tienen la misma secuencia idéntica, entonces se consideran duplicadas y se fusionarán.
En la salida de pantalla hay dos columnas:
- el primero es el número de secuencias caracterizadas y
- el segundo es el número de secuencias únicas restantes.
Entonces, después de ejecutar unique.seqs, pasamos de 125.540 en total a 20.469 secuencias únicas. Esto hará que nuestra vida sea mucho más fácil de seguir adelante.
Simplify names and groups with counts
Another thing to do to make our lives easier is to simplify the names and group files. mothur carefully keeps track of the name (i.e. sequence ID representing a set of uniques, .names file) and group (i.e. sample name, .groups file) for every single sequence. To make later commands shorter, we combine the data in these files into a single “count” file (.count_table).
So we’ll run count.seqs to generate a table where the rows are the names of the unique sequences and the columns are the names of the groups. The table is then filled with the number of times each unique sequence shows up in each group.
count.seqs
mothur > count.seqs(name=example.trim.contigs.good.names, group=example.contigs.good.groups)Using 2 processors.
It took 1 secs to create a table for 125540 sequences.
Total number of sequences: 125540
Output File Names:
/home/mcox/MiSeq_SOP/example.trim.contigs.good.count_tableThis will generate a file called example.trim.contigs.good.count_table. In subsequent commands, we’ll use this file with the count= option.
Let’s look at another summary of our data using this count table.
summary.seqs with count_table
mothur > summary.seqs(count=example.trim.contigs.good.count_table)Using 2 processors.
Start End NBases Ambigs Polymer NumSeqs
Minimum: 1 250 250 0 3 1
2.5%-tile: 1 252 252 0 3 3139
25%-tile: 1 253 253 0 3 31386
Median: 1 253 253 0 4 62771
75%-tile: 1 253 253 0 4 94156
97.5%-tile: 1 254 254 0 5 122402
Maximum: 1 294 294 0 8 125540
Mean: 1 252.889 252.889 0 3.82739
# of unique seqs: 20469
total # of seqs: 125540
Output File Names:
/home/mcox/MiSeq_SOP/example.trim.contigs.good.unique.summary
It took 0 secs to summarize 125540 sequences.We see that now, by including a fasta and count_table, we see “# of unique seqs” and “total # of seqs”. We would get the same result with summary.seqs(name=example.trim.contigs.good.names, group=example.contigs.good.groups).
Simplifica los nombres y grupos con los recuentos Otra cosa que debemos hacer para facilitar nuestras vidas es simplificar los archivos names ygroup. Motur realiza un seguimiento cuidadoso del nombre (* es decir, * ID de secuencia que representa un conjunto de archivos únicos, .names) y grupo (* es decir, * nombre de muestra, archivo.groups) para cada secuencia individual. Para acortar los comandos posteriores, combinamos los datos de estos archivos en un solo archivo de “conteo” (.count_table).
Entonces ejecutaremos count.seqs to para generar una tabla donde las filas son los nombres de las secuencias únicas y las columnas son los nombres de los grupos. La tabla se llena con la cantidad de veces que aparece cada secuencia única en cada grupo.
** count.seqs **
Esto generará un archivo llamado example.trim.contigs.good.count_table. En los comandos posteriores, usaremos este archivo con la opción count =.
Veamos otro resumen de nuestros datos usando esta tabla de recuento.
** summary.seqs concount_table **
Vemos que ahora, al incluir un fasta ** y ** count_table, vemos “# de seqs únicos” y “número total de seqs”. Obtendremos el mismo resultado con summary.seqs (name = example.trim.contigs.good.names, group = example.contigs.good.groups).
Alignment
Now we need to align our sequences to a reference alignment using align.seqs. Your sequences aren’t just random strings of nucleotides; they are part of a known gene (usually 16S or 18S). For ITS, SEE NOTE BELOW
We have two options for 16S alignment: Silva or GreenGenes. Only use Silva for alignment. GreenGenes contains incomplete (not full length) 16S sequences which make alignment difficult and error-prone. These databases are formatted specifically for mothur and can be found on the mothur wiki (Silva: https://www.mothur.org/wiki/Silva_reference_files, GreenGenes: https://www.mothur.org/wiki/Greengenes-formatted_databases).
We will also improve our alignment by including flip=T so that mothur checks both the forward and reverse complement of our sequences.
The first time you work with an amplicon (i.e. primer set), it’s a good idea to align to the full database to find out where your sequences align to. This is demonstrated below for V4 but should not be run during the workshop as it will take a long time and use a lot of RAM.
Instead, we will start aligning to the pre-made custom reduced alignment database you transferred earlier (silva.v4.align).
align.seqs
mothur > align.seqs(fasta=example.trim.contigs.good.unique.fasta, reference=silva.v128.v4.align, flip=T)Using 2 processors.
Reading in the /home/mcox/MiSeq_SOP/silva.v128.v4.align template sequences...
It took 35 to read 190661 sequences.
Aligning sequences from example.trim.contigs.good.unique.fasta ...
100
200
[truncated]
10200
10225
It took 66 secs to align 20469 sequences.
Output File Names:
/home/mcox/MiSeq_SOP/example.trim.contigs.good.unique.align
/home/mcox/MiSeq_SOP/example.trim.contigs.good.unique.align.reportNote: For fungal sequencing, it is best to skip the alignment step here. The ITS is too variable in length to generate a good alignment. Instead, we do an internal “needleman” alignment later at the precluster step. So if you are working with a fungal data set, skip to precluster at this point.
———————————-DO NOT RUN the following block———————————-
Alineación Ahora necesitamos alinear nuestras secuencias a una alineación de referencia usando align.seqs. Tus secuencias no son solo cadenas aleatorias de nucleótidos; son parte de un gen conocido (generalmente 16S o 18S). Para ITS, ** VER NOTA A CONTINUACIÓN **
Tenemos dos opciones para la alineación 16S: Silva o GreenGenes. ** Solo usa Silva para la alineación **. GreenGenes contiene secuencias 16S incompletas (no completas) que dificultan la alineación y son propensas a errores. Estas bases de datos están formateadas específicamente para mothur y se pueden encontrar en la wiki de mothur (Silva: https://www.mothur.org/wiki/Silva_reference_files, GreenGenes: https://www.mothur.org/wiki/Greengenes-formatted_databases) .
También mejoraremos nuestra alineación incluyendo flip = T para que mothur compruebe el complemento directo e inverso de nuestras secuencias.
La primera vez que trabaje con un amplicón (* i.e. * conjunto de cebadores), es una buena idea alinearse con la base de datos completa para saber dónde se alinean sus secuencias. Esto se demuestra a continuación para V4 pero ** no se debe ejecutar durante el taller **, ya que llevará mucho tiempo y usará mucha RAM.
En cambio, comenzaremos a alinearnos a la base de datos de alineación reducida personalizada prefabricada que transfirió anteriormente (silva.v4.align).
** align.seqs **
** Nota **: para la secuenciación de hongos, es mejor omitir el paso de alineación aquí. El ITS tiene una longitud demasiado variable para generar una buena alineación. En cambio, hacemos una alineación interna de “agujas” más adelante en el paso precluster. Entonces, si está trabajando con un conjunto de datos fúngicos, salte a precluster en este punto.
** ———————————- NO EJECUTE el siguiente bloque ——– ————————– **
Custom reduced database
Even with a reduced database, this may take several minutes to run. We will go over how this reduced database was created in the interim. First, we aligned to the full Silva database.
DO NOT RUN THIS
mothur > align.seqs(fasta=example.trim.contigs.good.unique.fasta, reference=silva.nr_v128.align, flip=T)Reading in the silva.nr_v128.align template sequences... DONE.
It took 506 to read 172418 sequences.
Aligning sequences from example.trim.contigs.good.unique.fasta ...
It took 363 secs to align 20469 sequences.
Output File Names:
example.trim.contigs.good.unique.align
example.trim.contigs.good.unique.align.reportThen we summary.seqs to see where we aligned to:
DO NOT RUN THIS
mothur > summary.seqs(fasta=example.trim.contigs.good.unique.align, count=example.trim.contigs.good.count_table) | Start | End | NBases | Ambigs | Polymer | NumSeqs
---------------- | ------- | ------- | ------- | -------- | ------- | -------
Minimum: | 13125 | 23439 | 250 | 0 | 3 | 1
2.5%-tile: | 13862 | 23444 | 252 | 0 | 3 | 3139
25%-tile: | 13862 | 23444 | 253 | 0 | 3 | 31386
Median: | 13862 | 23444 | 253 | 0 | 4 | 62771
75%-tile: | 13862 | 23444 | 253 | 0 | 4 | 94156
97.5%-tile: | 13862 | 23444 | 254 | 0 | 5 | 122402
Maximum: | 13875 | 26147 | 294 | 0 | 8 | 125540
Mean: | 13862 | 23444.1 | 252.889 | 0 | 3.82739 |
# of Seqs: | 20469
total # of seqs: | 125540So what does this mean?
You’ll see that the bulk of the sequences start at position 13862 (Start column) and end at position 23444 (End column). But some sequences start or end at other places. These deviants are likely due to an insertion or deletion at the terminal ends of the alignments. Sometimes you’ll see sequences that start and end at the same position indicating a very poor alignment, which is generally due to non-specific amplification.
It is these Start and End values that we need to make a custom reduced database. These values are different for different amplicons/primers and sometimes for different versions/releases of a database. So, once you know them for your amplicon and a given version of Silva, you won’t have to do this step again.
We make a database customized to our region of interest using the pcr.seqs command. To run this command, you need to have the full reference database (i.e. silva.nr_v128.align) and know where in that alignment your amplicon starts and ends (from the above summary). To remove the leading and trailing dots, we will set keepdots to false.
DO NOT RUN THIS
mothur > pcr.seqs(fasta=silva.nr_v128.align, start=13862, end=23444, keepdots=F)Output File Names:
silva.nr_v128.pcr.alignWe would then rename this file so we know which region we selected
DO NOT RUN THIS
mothur > system(mv silva.nr_v128.pcr.align silva.v128.v4.align)Note: Windows users should use “rename” rather than “mv”" since the system command calls upon the operating system of the computer. For this SOP, we are all running on Linux on Kabré.
Now we have a customized reduced reference alignment to align our sequences to. The nice thing about this reference is that instead of being 50,000 columns wide, it is now 13,425 columns wide which will save our hard drive some space and improve the overall alignment quality and efficiency.
silva.v4.align is the database we started an alignment to earlier. Once that is done, we summary.seqs again.
—————————————————————————–
Base de datos reducida personalizada Incluso con una base de datos reducida, esto puede tardar varios minutos en ejecutarse. Vamos a repasar cómo esta base de datos reducida se creó en el ínterin. Primero, nos alineamos con la base de datos completa de Silva.
** NO EJECUTE ESTO **
Luego hacemos summary.seqs para ver dónde nos alineamos:
** NO EJECUTE ESTO **
Entonces, ¿qué significa esto?
Verá que la mayor parte de las secuencias comienzan en la posición ** 13862 ** (columna de inicio) y terminan en la posición ** 23444 ** (columna final). Pero algunas secuencias comienzan o terminan en otros lugares. Estos desvíos probablemente se deben a una inserción o eliminación en los extremos terminales de las alineaciones. Algunas veces verá secuencias que comienzan y terminan en la misma posición, lo que indica una alineación muy pobre, que generalmente se debe a una amplificación no específica.
Son estos valores * Start * y * End * los que necesitamos para crear una base de datos reducida personalizada. Estos valores son diferentes para diferentes amplicones / cebadores y, en ocasiones, para diferentes versiones / lanzamientos de una base de datos. Entonces, una vez que los conozcas para tu amplicón y una versión dada de Silva, no tendrás que volver a hacer este paso.
Hacemos una base de datos personalizada para nuestra región de interés utilizando el comando pcr.seqs. Para ejecutar este comando, necesita tener la base de datos de referencia completa (* i.e. * silva.nr_v128.align) y saber dónde comienza y dónde termina su amplicón (del resumen anterior). Para eliminar los puntos inicial y final, estableceremos keepdots en falso.
** NO EJECUTE ESTO **
Entonces, cambiaríamos el nombre de este archivo para saber qué región seleccionamos
** Nota **: los usuarios de Windows deben usar “rename” en lugar de “mv” “ya que el comando del sistema llama al sistema operativo de la computadora. Para este SOP, todos estamos ejecutando Linux en Kabré.
Ahora tenemos ** una alineación de referencia reducida ** para alinear nuestras secuencias. Lo bueno de esta referencia es que en lugar de tener 50,000 columnas de ancho, ahora tiene 13,425 columnas de ancho que ahorrarán espacio en nuestro disco duro y mejorarán la calidad y eficiencia de la alineación general.
silva.v128.v4.align es la base de datos a la que comenzamos una alineación anterior. Una vez hecho esto, volvemos a summary.seqs.
** ———————————————— —————————– **
Post-alignment summary
summary.seqs
mothur > summary.seqs(fasta=example.trim.contigs.good.unique.align, count=example.trim.contigs.good.count_table)Using 2 processors.
Start End NBases Ambigs Polymer NumSeqs
Minimum: 8 9577 231 0 3 1
2.5%-tile: 8 9582 251 0 3 3139
25%-tile: 8 9582 252 0 3 31386
Median: 8 9582 252 0 4 62771
75%-tile: 8 9582 252 0 4 94156
97.5%-tile: 8 9582 253 0 5 122402
Maximum: 13 9582 256 0 8 125540
Mean: 8.00009 9582 251.886 0 3.82739
# of unique seqs: 20469
total # of seqs: 125540
Output File Names:
/home/mcox/MiSeq_SOP/example.trim.contigs.good.unique.summary
It took 2 secs to summarize 125540 sequences.
Resumen posterior a la alineación
** summary.seqs **
Rescreen
To make sure that everything overlaps the same region we’ll re-run screen.seqs. Note that we need to include the count table in this command so that we can update the table for the sequences we’re removing.
We’re also using the summary file we just made so we don’t have to figure out all the start and stop positions again. To this end, we will specify the following within the command: start=8, end=9582. These numbers are very different from those using the full Silva database because the numbering starts from the 5 prime end and we’ve removed a bunch of bases from that end with pcr.seqs
screen.seqs post alignment
mothur > screen.seqs(fasta=example.trim.contigs.good.unique.align, count=example.trim.contigs.good.count_table, summary=example.trim.contigs.good.unique.summary, start=8, end=9582)Using 2 processors.
Processing sequence: 100
Processing sequence: 200
[truncated]
Processing sequence: 10235
Processing sequence: 10234
Output File Names:
/home/mcox/MiSeq_SOP/example.trim.contigs.good.unique.good.summary
/home/mcox/MiSeq_SOP/example.trim.contigs.good.unique.good.align
/home/mcox/MiSeq_SOP/example.trim.contigs.good.unique.bad.accnos
/home/mcox/MiSeq_SOP/example.trim.contigs.good.good.count_table
It took 7 secs to screen 20469 sequences.And we summarize to see how re-screening impacted the data.
summary.seqs
mothur > summary.seqs(fasta=example.trim.contigs.good.unique.good.align, count=example.trim.contigs.good.good.count_table)Using 2 processors.
Start End NBases Ambigs Polymer NumSeqs
Minimum: 8 9582 231 0 3 1
2.5%-tile: 8 9582 251 0 3 3133
25%-tile: 8 9582 252 0 3 31322
Median: 8 9582 252 0 4 62644
75%-tile: 8 9582 252 0 4 93966
97.5%-tile: 8 9582 253 0 5 122155
Maximum: 8 9582 256 0 8 125287
Mean: 8 9582 251.888 0 3.82689
# of unique seqs: 20290
total # of seqs: 125287
Output File Names:
/home/mcox/MiSeq_SOP/example.trim.contigs.good.unique.good.summary
It took 2 secs to summarize 125287 sequences.
Rescreen Para asegurarnos de que todo se superpone a la misma región, volveremos a ejecutar screen.seqs. Tenga en cuenta que tenemos que incluir la tabla de recuento en este comando para que podamos actualizar la tabla de las secuencias que estamos eliminando.
También usamos el archivo de resumen que acabamos de crear para no tener que descifrar todas las posiciones de inicio y detención nuevamente. Para este fin, especificaremos lo siguiente dentro del comando: start = 8, end = 9582. Estos números son ** muy diferentes ** de los que usan la base de datos completa de Silva porque la numeración comienza desde el 5 prime y hemos eliminado un montón de bases de ese extremo con pcr.seqs
** post-alineación de screen.seqs **
Y resumimos para ver cómo la nueva evaluación impactó en los datos.
** summary.seqs **
Filter
Now that we know our sequences overlap the same alignment coordinates, we also want to make sure they only overlap that region. We’ll filter the sequences to remove the overhangs at both ends with vertical=T. Since we’ve done paired-end sequencing, this shouldn’t be much of an issue.
Also, when aligning, mothur puts in a lot of spaces and dots to show insertions and deletions compared to the alignment database. Therefore, there are many columns in the alignment that only contain gap characters (i.e. -). This is confusing to later commands (and causes them to fail) so we remove them all with trump=..
filter.seqs
mothur > filter.seqs(fasta=example.trim.contigs.good.unique.good.align, vertical=T, trump=.)Using 2 processors.
Creating Filter...
100
200
[truncated]
Running Filter...
100
200
[truncated]
2000
2028
Length of filtered alignment: 416
Number of columns removed: 9166
Length of the original alignment: 9582
Number of sequences used to construct filter: 20290
Output File Names:
example.filter
example.trim.contigs.good.unique.good.filter.fastaIF THIS FAILS
mothur > system(cp /home/capacitaciones/bioinformatica-UCR/example.trim.contigs.good.unique.good.filter.fasta ~/MiSeq_SOP)This means that our initial alignment was 9582 columns wide and that we were able to remove 9166 terminal gap characters using trump=. and vertical=T. The final alignment length is 416 columns.
Because we’ve perhaps created some redundancy across our sequences by trimming the ends, we can re-run unique.seqs
unique.seqs again
mothur > unique.seqs(fasta=example.trim.contigs.good.unique.good.filter.fasta, count=example.trim.contigs.good.good.count_table)1000 997
2000 1992
[truncated]
20000 19922
20290 20212
Output File Names:
example.trim.contigs.good.unique.good.filter.count_table
example.trim.contigs.good.unique.good.filter.unique.fastaThis identified 78 duplicate sequences that we’ve now merged with previous unique sequences.
Filtrar Ahora que sabemos que nuestras secuencias se superponen con las mismas coordenadas de alineación, también queremos asegurarnos de que solo se superpongan a esa región. Filtraremos las secuencias para eliminar los voladizos en ambos extremos con vertical = T. Como hemos hecho la secuencia de emparejamiento final, esto no debería ser un gran problema.
Además, al alinear, mothur coloca muchos espacios y puntos para mostrar las inserciones y eliminaciones en comparación con la base de datos de alineación. Por lo tanto, hay muchas columnas en la alineación que solo contienen caracteres de espacio (* i.e. * -). Esto es confuso para los comandos posteriores (y hace que fallen) por lo que los eliminamos todos con trump = ..
** filter.seqs **
Esto significa que nuestra alineación inicial fue de 9582 columnas de ancho y que pudimos eliminar 9166 caracteres de brecha terminal usando trump = . yvertical = T. La longitud de alineación final es de 416 columnas.
Debido a que tal vez hayamos creado alguna redundancia en nuestras secuencias al recortar los extremos, podemos volver a ejecutar unique.seqs
** unique.seqs de nuevo **
Esto identificó 78 secuencias duplicadas que ahora hemos fusionado con secuencias únicas anteriores.
Pre-cluster
The next thing we want to do to further de-noise our sequences is to pre-cluster the sequences that are very similar. We’ll use the pre.cluster command allowing for up to 2 differences between sequences by adding diffs=2 to the command.
This command will split the sequences by group and then sort them by abundance and go from most abundant to least and identify sequences that are within 2 nt of each other. If they are, then they get merged.
We generally favor allowing 1 difference for every 100 bp of sequence, hence a 250bp amplicon is diffs=2:
pre.cluster
mothur > pre.cluster(fasta=example.trim.contigs.good.unique.good.filter.unique.fasta, count=example.trim.contigs.good.unique.good.filter.count_table, diffs=2)Using 2 processors.
Processing group L62:
0 2625 34
Total number of sequences before pre.cluster was 2659.
pre.cluster removed 1349 sequences.
It took 1 secs to cluster 2659 sequences.
Processing group S62:
100 1991 668
200 1825 834
[truncated]
18700 5002 13737
18739 5001 13738
Total number of sequences before pre.cluster was 18739.
pre.cluster removed 13738 sequences.
It took 5 secs to cluster 18739 sequences.
It took 6 secs to run pre.cluster.
Output File Names:
example.trim.contigs.good.unique.good.filter.unique.precluster.fasta
example.trim.contigs.good.unique.good.filter.unique.precluster.count_table
example.trim.contigs.good.unique.good.filter.unique.precluster.L62.map
example.trim.contigs.good.unique.good.filter.unique.precluster.S62.mapNote: For fungal data, you should include “align=needleman” here in the precluster command to do the internal alignment
summary.seqs
mothur > summary.seqs(fasta=example.trim.contigs.good.unique.good.filter.unique.precluster.fasta, count=example.trim.contigs.good.unique.good.filter.unique.precluster.count_table)Using 2 processors.
Start End NBases Ambigs Polymer NumSeqs
Minimum: 1 416 231 0 3 1
2.5%-tile: 1 416 251 0 3 3133
25%-tile: 1 416 252 0 3 31322
Median: 1 416 252 0 4 62644
75%-tile: 1 416 252 0 4 93966
97.5%-tile: 1 416 253 0 5 122155
Maximum: 1 416 256 0 8 125287
Mean: 1 416 251.896 0 3.82339
# of unique seqs: 5619
total # of seqs: 125287
Output File Names:
example.trim.contigs.good.unique.good.filter.unique.precluster.summary
It took 0 secs to summarize 125287 sequences.Note that we’ve gone from 20,208 to 5,614 unique sequences but our total # of sequences has not changed. This is because no sequences were removed; they were just grouped similarly to unique.seqs but allowing up to 2 differences.
It is not uncommon to drastically reduce your unique sequence count with this command. The MiSeq is good but not perfect.
At this point we have removed as much sequencing error as we can, and it is time to turn our attention to removing chimeras.
Pre-cluster Lo siguiente que queremos hacer para seguir reduciendo el ruido de nuestras secuencias es agrupar previamente las secuencias que son muy similares. Usaremos el comando pre.cluster que permite ** hasta 2 diferencias ** entre secuencias al agregardiffs = 2 al comando.
Este comando dividirá las secuencias por grupo y luego las clasificará por abundancia e irá de más abundante a menos e identificará las secuencias que están dentro de 2 nt una de la otra. Si lo son, se fusionan.
En general, estamos a favor de permitir 1 diferencia por cada 100 pb de secuencia, por lo tanto, un amplicón de 250 pb es diffs = 2:
** pre.cluster **
** Nota: ** Para datos de hongos, debe incluir “align = needleman” aquí en el comando precluster para hacer la alineación interna
** summary.seqs **
Tenga en cuenta que hemos pasado de 20,208 a 5,614 secuencias únicas, pero nuestro número total de secuencias no ha cambiado. Esto se debe a que no se eliminaron las secuencias; simplemente se agruparon de manera similar a unique.seqs pero permitiendo hasta 2 diferencias.
No es raro reducir drásticamente el recuento de secuencias único con este comando. El MiSeq es bueno pero no perfecto.
En este punto, hemos eliminado todo el error de secuenciación que hemos podido, y es hora de centrar nuestra atención en eliminar las quimeras.
Chimeras
Chimeras occur when the polymerase falls off one sequence and then reattaches to another during a single round of PCR. This results in a sequence that looks unique but is, in fact, just half of one and half of another e.g. not real.
We will remove chimeras with the UCHIME algorithm that is called within mothur using the chimera.uchime command.
Again, this command will split the data by sample and check for chimeras.
Our preferred way of doing this is to use the abundant sequences as our reference. In addition, if a sequence is flagged as chimeric in one sample, the default (dereplicate=F) is to remove it from all samples. Our experience suggests that this is a bit aggressive since we’ve seen rare sequences get flagged as chimeric when they’re the most abundant sequence in another sample.
chimera.uchime
mothur > chimera.uchime(fasta=example.trim.contigs.good.unique.good.filter.unique.precluster.fasta, count=example.trim.contigs.good.unique.good.filter.unique.precluster.count_table, dereplicate=t)uchime by Robert C. Edgar
http://drive5.com/uchime
This code is donated to the public domain.
Checking sequences from example.trim.contigs.good.unique.good.filter.unique.precluster.fasta ...
uchime v4.2.40
uchime v4.2.40
by Robert C. Edgar
by Robert C. Edgar
http://drive5.com/uchime
http://drive5.com/uchime
This code is donated to the public domain.
00:00 18Mb 0.1% 0Readi0ng ex:0am0 pl 1e.8Mtrb i m. 0co.1nt% igRs.eagodiodng.u niexquame.plgoeod.fil.tterir.m.uncoique.nprtiec00:00 19Mb 100.0% Reading example.trim.contigs.good.unique.good.filter.unique.precluster.temp
00:00 19Mb 1310 sequences
00:00 20Mb 100.0% Reading example.trim.contigs.good.unique.good.filter.unique.precluster.temp27375.temp
00:00 20Mb 5001 sequences
00:09 2.9Mb 100.0% 78/1309 chimeras found (6.0%) nd (3. 5% 9)
It took 9 secs to check 1310 sequences from group L62.
01:39 5.2Mb 100.0% 1266/5000 chimeras found (25.3%)
It took 100 secs to check 5001 sequences from group S62.
Output File Names:
example.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.count_table
example.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.chimeras
example.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.accnosThis command only identifies chimeras. We then need to remove them from the data using remove.seqs
remove.seqs
mothur > remove.seqs(fasta=example.trim.contigs.good.unique.good.filter.unique.precluster.fasta, count=example.trim.contigs.good.unique.good.filter.unique.precluster.count_table, accnos=example.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.accnos)[NOTE]: The count file should contain only unique names, so mothur assumes your fasta, list and taxonomy files also contain only uniques.
Removed 1297 sequences from your fasta file.
Removed 3477 sequences from your count file.
Output File Names:
example.trim.contigs.good.unique.good.filter.unique.precluster.pick.fasta
example.trim.contigs.good.unique.good.filter.unique.precluster.pick.count_tablesummary.seqs
mothur > summary.seqs(fasta=example.trim.contigs.good.unique.good.filter.unique.precluster.pick.fasta, count=example.trim.contigs.good.unique.good.filter.unique.precluster.pick.count_table)Using 2 processors.
Start End NBases Ambigs Polymer NumSeqs
Minimum: 1 416 241 0 3 1
2.5%-tile: 1 416 251 0 3 3046
25%-tile: 1 416 252 0 3 30453
Median: 1 416 252 0 4 60906
75%-tile: 1 416 252 0 4 91358
97.5%-tile: 1 416 253 0 5 118765
Maximum: 1 416 255 0 8 121810
Mean: 1 416 251.894 0 3.81972
# of unique seqs: 4322
total # of seqs: 121810
Output File Names:
example.trim.contigs.good.unique.good.filter.unique.precluster.pick.summary
It took 0 secs to summarize 121810 sequences.Note that we went from 125,287 to 121,810 sequences for a reduction of about 3%; this is a reasonable number of sequences to be flagged as chimeric. In the end, the percent chimeras is a result of the polymerase you used in PCR. A high-fidelity polymerase is less likely to fall off and create chimeras and fewer cycles also result in fewer chimeras.
Quimeras Las quimeras se producen cuando la polimerasa se desprende de una secuencia y luego se vuelve a unir a otra durante una única ronda de PCR. Esto da como resultado una secuencia que se ve única pero que, de hecho, es solo la mitad de una y la mitad de otra * por ej. * No es real.
Eliminaremos quimeras con el algoritmo * UCHIME * que se llama dentro de mothur utilizando el comando chimera.uchime.
De nuevo, este comando dividirá los datos por muestra y comprobará quimeras.
Nuestra forma preferida de hacer esto es usar las abundantes secuencias como nuestra referencia. Además, si una secuencia se marca como quimérica en una muestra, el valor predeterminado (dereplicate = F) es eliminarla de todas las muestras. Nuestra experiencia sugiere que esto es un poco agresivo ya que hemos visto secuencias raras que se marcan como quiméricas cuando son la secuencia más abundante en otra muestra.
** quimera.uchime **
Este comando solo identifica quimeras. Entonces necesitamos eliminarlos de los datos usando remove.seqs
** remove.seqs **
** summary.seqs **
Tenga en cuenta que pasamos de 125,287 a 121,810 secuencias para una reducción de aproximadamente 3%; este es un número razonable de secuencias que se marcarán como quiméricas. Al final, el porcentaje de quimeras es el resultado de la polimerasa que usaste en la PCR. Una polimerasa de alta fidelidad tiene menos probabilidades de caer y crear quimeras y menos ciclos también dan como resultado menos quimeras.
Undesirables
As a final quality control step, we need to see if there are any “undesirables” in our data set. Even after all of this cleanup, you can still have “good” sequences that are actually contaminants.
Sometimes when we pick a primer set, they will amplify other stuff that gets to this point in the pipeline such as 18S rRNA gene fragments or 16S rRNA from Archaea, chloroplasts, and mitochondria.
There’s also just the random stuff that we want to get rid of. Now you may say, “But wait I want that stuff”. Fine. But, the primers we use, are only supposed to amplify members of the Bacteria and if they’re hitting Eukaryota or Archaea, then its a mistake. It is important to note that if you do keep these sequences, they do not represent the overall archaeal or eukaryotic community. It would be incorrect to use this data as a shortcut to not do an archaea or eukaryote specific MiSeq run.
You may also have other “undesirables” specific to your system of study. For example, the rumen is not home to any Cyanobacteria as it is a dark, anaerobic environment. However, we always get Cyanobacteria sequences back because unfortunately, chloroplasts (very abundant in the feed of herbivores!) look a lot like Cyanobacteria. So we remove them. Not every researcher does this but it is an accepted step as long as you note it in your methods section.
The first step is to classify your sequences so you know what to remove. Let’s go ahead and classify using the Bayesian classifier with the classify.seqs command. You can classify to either Silva or GreenGenes at this point and it will make very little difference. We will use Silva herefor bacterial sequences. We will also specific a cutoff of 80. This roughly means that we are 80%+ confident in our classification, a good thing if we’re going to remove groups. We want to be sure of what we’re removing.
classify.seqs
mothur > classify.seqs(fasta=example.trim.contigs.good.unique.good.filter.unique.precluster.pick.fasta, count=example.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.count_table, reference=silva.v128.v4.align, taxonomy=silva.nr_v128.tax, cutoff=80)Using 2 processors.
Generating search database... DONE.
It took 38 seconds generate search database.
Reading in the silva.nr_v128.tax taxonomy... DONE.
Calculating template taxonomy tree... DONE.
Calculating template probabilities... DONE.
Classifying sequences from example.trim.contigs.good.unique.good.filter.unique.precluster.pick.fasta ...
Processing sequence: 100
Processing sequence: 100
[truncated]
[WARNING]: M03662_8_000000000-AK09K_1_2104_20318_19579 could not be classified. You can use the remove.lineage command with taxon=unknown; to remove such sequences.
[truncated]
Processing sequence: 2164
Processing sequence: 2158
It took 62 secs to classify 4322 sequences.
It took 1 secs to create the summary file for 4322 sequences.
Output File Names:
example.trim.contigs.good.unique.good.filter.unique.precluster.pick.nr_v128.wang.taxonomy
example.trim.contigs.good.unique.good.filter.unique.precluster.pick.nr_v128.wang.tax.summaryNote: Silva and GreenGenes are the options provided, in correct format for use in alignment and classification in mothur, by Pat Schloss and the mothur developers. GreenGenes is appropriate for most well-characterized host-associated gut communities, but is poor for environmental data sets. If you are working with environmental data, it maybe a good idea to find or make an alternative database in order to get better classifications. See an example for lake systems created by the McMahon lab here on campus: https://github.com/McMahonLab/TaxAss/tree/master/FreshTrain-files.
Note: With most data sets, you will get warnings about things that could not be classified. These are labeled as domain ‘unknown’ if the classifier cannot classify the sequence to Bacteria, Archaea, or Eukaryota.
Now that everything is classified, we want to remove our undesirables. We do this with the remove.lineage command. We’ll just remove non-Bacteria domains and completely unclassifieds (“unknown”).
remove.lineage
mothur > remove.lineage(fasta=example.trim.contigs.good.unique.good.filter.unique.precluster.pick.fasta, count=example.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.count_table, taxonomy=example.trim.contigs.good.unique.good.filter.unique.precluster.pick.nr_v128.wang.taxonomy, taxon=unknown;-Archaea;-Eukaryota;)[NOTE]: The count file should contain only unique names, so mothur assumes your fasta, list and taxonomy files also contain only uniques.
Output File Names:
example.trim.contigs.good.unique.good.filter.unique.precluster.pick.nr_v128.wang.pick.taxonomy
example.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.fasta
example.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.pick.count_tableHere is an example including the removal of full lineages (Cyanobacteria and mitochondria), as we would do in the rumen. NOTE that the hyphens are list separators, and the taxonomic levels are separated by semicolons.
mothur > remove.lineage(fasta=example.trim.contigs.good.unique.good.filter.unique.precluster.pick.fasta, count=example.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.count_table, taxonomy=example.trim.contigs.good.unique.good.filter.unique.precluster.pick.nr_v128.wang.taxonomy, taxon=unknown;-Archaea;-Eukaryota;-Bacteria;Cyanobacteria;-Bacteria;Proteobacteria;Alphaproteobacteria;Rickettsiales;mitochondria;) summary.seqs
mothur > summary.seqs(fasta=example.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.fasta, count=example.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.pick.count_table)Using 2 processors.
Start End NBases Ambigs Polymer NumSeqs
Minimum: 1 416 241 0 3 1
2.5%-tile: 1 416 251 0 3 3015
25%-tile: 1 416 252 0 3 30143
Median: 1 416 252 0 4 60285
75%-tile: 1 416 252 0 4 90427
97.5%-tile: 1 416 253 0 5 117555
Maximum: 1 416 255 0 8 120569
Mean: 1 416 251.892 0 3.81681
# of unique seqs: 4303
total # of seqs: 120569
Output File Names:
example.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.summary
It took 0 secs to summarize 120569 sequences.At this point, we have cleaned-up our data as far as possible. Our file names are very long at this point, so let’s copy and rename them before moving onto some analysis.
Rename final files
mothur > system(cp example.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.fasta example.final.fasta)
mothur > system(cp example.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.pick.count_table example.final.count_table)
Indeseable Como paso final de control de calidad, necesitamos ver si hay algún “indeseable” en nuestro conjunto de datos. Incluso después de toda esta limpieza, aún puede tener secuencias “buenas” que en realidad son contaminantes.
A veces, cuando elegimos un conjunto de cebadores, amplifican otras cosas que llegan a este punto en la tubería, como los fragmentos del gen 18S rRNA o 16S rRNA de Archaea, cloroplastos y mitocondrias.
También hay cosas al azar de las que queremos deshacernos. Ahora puedes decir: “Pero espera, quiero esas cosas”. Multa. Pero, los cebadores que usamos, solo se supone que amplifican a los miembros de las Bacterias y si están golpeando Eukaryota o Archaea, entonces es un error. Es importante tener en cuenta que si mantiene estas secuencias, ** no representan ** a la comunidad global de arqueas o eucariotas. Sería incorrecto usar esta información como un atajo para no hacer una ejecución MiSeq específica de arqueas o eucariotas.
También puede tener otros “indeseables” específicos para su sistema de estudio. Por ejemplo, el rumen no alberga ninguna cianobacteria, ya que es un ambiente oscuro y anaeróbico. Sin embargo, siempre recuperamos las secuencias de las cianobacterias porque, desafortunadamente, los cloroplastos (¡muy abundantes en la alimentación de los herbívoros!) Se parecen mucho a las cianobacterias. Entonces los eliminamos. No todos los investigadores hacen esto, pero es un paso aceptado siempre que lo notes en la sección de métodos.
El primer paso es clasificar tus secuencias para que sepas qué eliminar. Avancemos y clasifiquemos usando el ** clasificador bayesiano ** con el comando classify.seqs. Puede clasificar a Silva o GreenGenes en este punto y hará muy poca diferencia. Utilizaremos * Silva * aquí para las secuencias bacterianas. También especificaremos un * punto de corte de 80 *. Esto significa aproximadamente que tenemos un 80% de confianza en nuestra clasificación, algo bueno si vamos a eliminar grupos. Queremos estar seguros de lo que estamos eliminando.
** Nota: ** Silva y GreenGenes son las opciones provistas, en el formato correcto para uso en alineación y clasificación en mothur, por Pat Schloss y los desarrolladores de mothur. GreenGenes es apropiado para la mayoría de las comunidades intestinales asociadas al host bien caracterizadas, pero es ** pobre ** para los conjuntos de datos ambientales. Si está trabajando con datos ambientales, tal vez sea una buena idea buscar o crear una base de datos alternativa para obtener mejores clasificaciones. Vea un ejemplo de sistemas de lagos creados por el laboratorio de McMahon aquí en el campus: https://github.com/McMahonLab/TaxAss/tree/master/FreshTrain-files.
** Nota: ** Con la mayoría de los conjuntos de datos, recibirá advertencias sobre cosas que no se pudieron clasificar. Estos se etiquetan como dominio “desconocido” si el clasificador no puede clasificar la secuencia en Bacteria, Archaea o Eukaryota.
Ahora que todo está clasificado, queremos eliminar nuestros indeseables. Hacemos esto con el comando remove.lineage. Simplemente eliminaremos los dominios que no sean de Bacteria y los completamente no clasificados (“desconocido”).
** remove.lineage **
** summary.seqs **
En este punto, hemos limpiado nuestros datos en la medida de lo posible. Nuestros nombres de archivo son muy largos en este punto, así que copiemos y cambiemos el nombre antes de pasar a un análisis.
** Renombrar archivos finales **
Define OTUs
Now we can define operational taxonomic units (OTUs), the sequencing proxy for a microbial species. We do this by calculating distances between sequences (how difference they are from each other) with dist.seqs and then clustering these distances based on a difference cutoff with cluster.split. In general, cutoffs are
- 0.03 = 3% different or 97% similar ~ species
- 0.05 = 5% different or 95% similar ~ genus
- 0.1 = 10% different or 90% similar ~ family
According to this paper.
Definir OTUs Ahora podemos definir unidades taxonómicas operativas (OTU), el proxy de secuencia para una especie microbiana. Hacemos esto calculando las distancias entre las secuencias (cómo se diferencian entre sí) con dist.seqs y luego agrupando estas distancias en función de un corte de diferencia concluster.split. En general, los puntos de corte son
- 0.03 = 3% especies diferentes o 97% similares ~
- 0.05 = 5% de género diferente o 95% similar ~
- 0.1 = 10% diferente o 90% similar ~ familia
Generate distance matrix
dist.seqs
mothur > dist.seqs(fasta=example.final.fasta)Using 2 processors.
0 0
100 0
[truncated]
4300 25
4302 25
Output File Names:
example.final.dist
It took 40 seconds to calculate the distances for 4303 sequences.Note: For fungal data, you should use pairwise.seqs instead of dist.seqs because your sequences are not aligned by align.seqs which is required by dist.seqs.
Note: There are a couple of tricks to decrease your distance matrix size for large data sets (I’ve seen V3-4 data get up to 1TB). You can set a “cutoff” in dist.seqs. This causes mothur to not save distance greater than the cutoff. For example, if you’re only going to look at genus and species level OTUs, you could set your cutoff to 0.1. You can also output in lower triangle form (output=lt), which only saves one distance for every pairwise comparison.
Generar matriz de distancia
** dist.seqs **
** Nota : para datos de hongos **, debe usar pairwise.seqs en lugar dedist.seqs porque sus secuencias no están alineadas por align.seqs que es requerido pordist.seqs .
** Nota **: Hay un par de trucos para disminuir el tamaño de la matriz de distancia para grandes conjuntos de datos (he visto datos de V3-4 que llegan a 1TB). Puede establecer un “punto de corte” en dist.seqs. Esto hace que mothur no guarde una distancia mayor que el límite. Por ejemplo, si solo va a mirar las OTU genus y especies, podría establecer su límite en 0.1. También puede dar salida en forma de triángulo inferior (salida = lt), que solo guarda una distancia por cada comparación por pares.
Cluster
For clustering, we have a number of options
- nearest neighbor
- furthest neighbor
- average neighbor
- opticlust
You can find full definitions of these methods here but the general view is that furthest is too strict (creates too many OTUs), nearest is too lenient (too few OTUs), and average is in the middle. Sort of a Goldilocks situation.
opticlust is a very new method developed by the Schloss group, and it appears to perform as good or better than average neighbor. For more information, see http://msphere.asm.org/content/2/2/e00073-17
We will use average here.
cluster.split
mothur > cluster.split(column=example.final.dist, count=example.final.count_table, method=average)Using 2 processors.
Splitting the file...
It took 211 seconds to split the distance file.
example.final.dist.0.temp
example.final.dist.31.temp
example.final.dist.18.temp
example.final.dist.45.temp
example.final.dist.19.temp
example.final.dist.34.temp
example.final.dist.38.temp
example.final.dist.22.temp
example.final.dist.37.temp
example.final.dist.36.temp
example.final.dist.56.temp
example.final.dist.27.temp
example.final.dist.47.temp
example.final.dist.58.temp
Reading example.final.dist.0.temp
Reading example.final.dist.41.temp
********************#****#****#****#****#****#****#****#****#****#****#
Reading matrix: ||||||||||||||||||||||||||||||||||||||||||||||||||||
***********************************************************************
Clustering example.final.dist.41.temp
[truncated]
Cutoff was 0.15 changed cutoff to 0.07
It took 2 seconds to cluster
Merging the clustered files...
It took 3 seconds to merge.
Output File Names:
example.final.an.unique_list.list
Racimo
Para la agrupación, tenemos una serie de opciones
- Vecino más cercano
- el vecino más lejano
- vecino promedio
- opticlust
Puede encontrar definiciones completas de estos métodos [aquí] (https://www.mothur.org/wiki/Cluster.split) pero la opinión general es que la más lejana es demasiado estricta (crea demasiadas OTU), la más cercana es demasiado indulgente ( muy pocas OTU), y el promedio está en el medio. Una especie de situación Goldilocks.
opticlust es un método muy nuevo desarrollado por el grupo Schloss, y parece funcionar tan bien o mejor que el vecino promedio. Para obtener más información, consulte http://msphere.asm.org/content/2/2/e00073-17
Usaremos el promedio aquí.
** cluster.split **
Generate OTU matrix
Now that we have OTUs, we’d like to look at them, right? Well, distance and cluster files are not for human eyes. So we use make.shared to combine these data with our sample names to create a table we can understand. We will have mothur only give us the species-level OTUs label=0.03 but you could ask for any level that you like (as long as it’s below the cutoff you may have used in dist.seqs).
make.shared
mothur > make.shared(list=example.final.an.unique_list.list, count=example.final.count_table, label=0.03)0.03
Output File Names:
example.final.an.unique_list.sharedLet’s look at it!!! This file has many, many columns and is difficult to look at in the Terminal. So, we will open the .shared file in Excel or a text editor. You’ll see that for each sample, we have how many times each OTU occurs in that sample.
| label | Group | numOtus | Otu0001 | Otu0002 … |
|---|---|---|---|---|
| 0.03 | L62 | 2204 | 189 | 678 … |
| 0.03 | S62 | 2204 | 7351 | 6438 … |
Generar matriz OTU
Ahora que tenemos OTU, nos gustaría verlos, ¿verdad? Bueno, los archivos de distancia y clúster no son para ojos humanos. Entonces usamos make.shared para combinar estos datos con nuestros nombres de muestra para crear una tabla que podamos entender. Tendremos que mothur solo nos proporcione las OTU de nivel de especie label = 0.03 pero puede solicitar cualquier nivel que desee (siempre que esté por debajo del límite que haya utilizado endist.seqs).
** make.shared **
¡¡Miremos eso !!! Este archivo tiene muchas, muchas columnas y es difícil de ver en la Terminal. Entonces, abriremos el archivo .shared en Excel o un editor de texto. Verá que para cada muestra, tenemos cuántas veces ocurre cada OTU en esa muestra.
| label | Group | numOtus | Otu0001 | Otu0002 … |
|---|---|---|---|---|
| 0.03 | L62 | 2204 | 189 | 678 … |
| 0.03 | S62 | 2204 | 7351 | 6438 … |
Classify OTUs
Now that you have all these OTUs, what are they? We find out by repeating classify.seqs on our final files and then grouping based on our OTUs with classify.otu. Here, we will use GreenGenes instead of Silva. Like classifying before, you can use either. In general, GreenGenes will give fewer unclassifieds and more classifications down to lower levels. But when you go to research those classifications, they could be from a clone library or other non-culturing technique so there is little you can say about that species. Silva is only full-length 16S so you have fewer sequences to compare to and get more unclassifieds. However, the classifications are more often to cultured species.
There are also specific taxonomic groups that are more or less represented in each database. This would be specific to your system and you’d have to classify to both to find out!
classify.seqs
mothur > classify.seqs(fasta=example.final.fasta, count=example.final.count_table, template=gg_13_8_99.fasta, taxonomy=gg_13_8_99.gg.tax, cutoff=80)Using 2 processors.
Generating search database... DONE.
It took 103 seconds generate search database.
Reading in the gg_13_8_99.gg.tax taxonomy... DONE.
Calculating template taxonomy tree... DONE.
Calculating template probabilities... DONE.
It took 283 seconds get probabilities.
Classifying sequences from example.final.fasta ...
Processing sequence: 100
Processing sequence: 100
[truncated]
Processing sequence: 2155
Processing sequence: 2148
It took 49 secs to classify 4303 sequences.
It took 1 secs to create the summary file for 4303 sequences.
Output File Names:
example.final.gg.wang.taxonomy
example.final.gg.wang.tax.summaryclassify.otu
mothur > classify.otu(list=example.final.an.unique_list.list, taxonomy=example.final.gg.wang.taxonomy, count=example.final.count_table, label=0.03, cutoff=80, basis=otu, probs=F)0.03 1812
Output File Names:
example.final.an.unique_list.0.03.cons.taxonomy
example.final.an.unique_list.0.03.cons.tax.summaryLooking at it.
mothur > system(head -5 example.final.an.unique_list.0.03.cons.taxonomy)OTU Size Taxonomy
Otu0001 21662 k__Bacteria;p__Bacteroidetes;c__Bacteroidia;o__Bacteroidales;f__Prevotellaceae;g__Prevotella;s__ruminicola;
Otu0002 7553 k__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Aeromonadales;f__Succinivibrionaceae;f__Succinivibrionaceae_unclassified;f__Succinivibrionaceae_unclassified;
Otu0003 5599 k__Bacteria;p__Firmicutes;c__Clostridia;o__Clostridiales;f__Veillonellaceae;g__Succiniclasticum;g__Succiniclasticum_unclassified;
Otu0004 3131 k__Bacteria;p__Proteobacteria;p__Proteobacteria_unclassified;p__Proteobacteria_unclassified;p__Proteobacteria_unclassified;p__Proteobacteria_unclassified;p__Proteobacteria_unclassified;Now, the species level classifications are dicey. Because we don’t have full length 16S, we can’t say with certainty that OTU2 is Prevotella ruminicola. What we can say is that OTU2 is most closely related to Prevotella ruminicola, so if we were going to guess at OTU2’s function in the cow gut, data on Prevotella ruminicola would be a good place to start. And that’s not even getting into strain level diversity.
Coverage
How to you know you have enough sequences to see all/most of the diversity in your samples? Do you need to redo some samples? You can visualize this by rarefaction and calculate it with Good’s coverage.
Classify OTUs Ahora que tienes todas estas OTU, ¿qué son? Descubrimos repitiendo classify.seqs en nuestros archivos finales y luego agrupamos en base a nuestras OTUs conclassify.otu. Aquí, usaremos ** GreenGenes ** en lugar de Silva. Al igual que clasificar antes, puede usar cualquiera. En general, GreenGenes otorgará menos clasificaciones no clasificadas y más a niveles inferiores. Pero cuando vayas a investigar esas clasificaciones, podrían ser de una biblioteca de clones u otra técnica que no sea de cultivo, por lo que es poco lo que puedes decir sobre esa especie. Silva solo tiene 16S de largo completo, por lo que tienes menos secuencias para comparar y obtener más clasificaciones sin clasificar. Sin embargo, las clasificaciones son más a menudo especies cultivadas.
También hay grupos taxonómicos específicos que están más o menos representados en cada base de datos. ¡Esto sería específico para su sistema y usted tendría que clasificar a ambos para averiguarlo!
** classify.seqs **
** classify.otu **
Mirándolo
Ahora, las clasificaciones del nivel de especie son inciertas. Debido a que no tenemos 16S de longitud completa, no podemos decir con certeza que OTU2 es * Prevotella ruminicola . Lo que podemos decir es que OTU2 está más estrechamente relacionado con Prevotella ruminicola , así que si íbamos a adivinar la función de OTU2 en el intestino de la vaca, los datos sobre Prevotella ruminicola * serían un buen lugar para comenzar. Y eso ni siquiera entra en la diversidad del nivel de tensión.
Cobertura ¿Cómo sabes que tienes suficientes secuencias para ver toda / la mayoría de la diversidad en tus muestras? ¿Necesitas volver a hacer algunas muestras? Puede visualizar esto por rarefacción y calcularlo con la cobertura de Good.
Rarefaction
rarefaction.single
mothur > rarefaction.single(shared=example.final.an.unique_list.shared)Using 2 processors.
Processing group L62
0.03
Processing group S62
0.03
Output File Names:
example.final.an.unique_list.groups.rarefactionThis will produce a file called example.final.an.unique_list.groups.rarefaction. Looking at it…
mothur > system(head -15 example.final.an.unique_list.groups.rarefaction)numsampled 0.03-L62 lci-L62 hci-L62 0.03-S62 lci-S62 hci-S62
1 1 1 1 1 1 1
100 71.4920 63 80 65.6890 57 74
200 123.8100 112 136 111.8440 101 123
300 167.9220 154 183 150.7890 137 165
400 206.6730 192 223 185.2670 169 200
500 241.3460 225 258 216.2650 198 232
600 272.5670 255 290 244.5300 225 262
700 300.9320 283 320 270.6930 251 290
800 327.0350 309 346 295.0480 275 316
900 351.3850 332 373 317.8860 297 340
1000 373.9970 354 395 339.1980 317 361
1100 395.0880 376 416 359.2910 337 382
1200 414.7530 395 436 378.2990 355 402
1300 433.4830 414 455 396.3390 373 420it’s difficult to tell what’s going on. The best way to use this data is to graph it and look from the new unique sequences gained per ## of sequences added to level off. At a certain point, you need 1000s of sequences to get only 1 more unique sequence so it’s not worth it (and probably not even real)
For example (made in Excel) 
However, this doesn’t numerically/statistically show you your coverage and is subject to axis-skewing. For example, our coverage looks pretty good in the above plot, but if we compress the x-axis, it doesn’t look nearly as good. 
Thus, it is better to also calculate coverage with another method, like Good’s.
Rarefacción
** rarefaction.single **
Esto producirá un archivo llamado example.final.an.unique_list.groups.rarefaction. Mirándolo …
es difícil saber qué está pasando. La mejor manera de utilizar estos datos es graficarlo y observar las nuevas secuencias únicas obtenidas por ## de secuencias añadidas para nivelar. En cierto punto, necesitas miles de secuencias para obtener solo una secuencia única más por lo que no vale la pena (y probablemente ni siquiera sea real)
Por ejemplo (hecho en Excel) ! [] (Rarefaction.png)
Sin embargo, esto no muestra numérica / estadísticamente su cobertura y está sujeto a sesgos de eje. Por ejemplo, nuestra cobertura se ve bastante bien en el gráfico anterior, pero si comprimimos el eje x, no parece tan bueno. ! [] (Rarefaction2.png)
Por lo tanto, es mejor también calcular la cobertura con otro método, como el de Good.
Good’s coverage
The metric is an estimate of the percent of diversity (species) captured in a sample by your current sequencing data. In general, we’re looking for 95% or better. We can also look at the number of sequence nseqs and number of unique OTUs sobs in each sample.
summary.single
mothur > summary.single(shared=example.final.an.unique_list.shared, label=0.03, calc=nseqs-sobs-coverage)Processing group L62
0.03
Processing group S62
0.03
Output File Names:
example.final.an.unique_list.groups.summaryIf we look at it,
mothur > system(head -3 example.final.an.unique_list.groups.summary)label group nseqs sobs coverage
0.03 L62 12597.000000 801.000000 0.993887
0.03 S62 107972.000000 1710.000000 0.997046we see that we have very good coverage of these samples, >99%
Cobertura “Good’s” La métrica es una estimación del porcentaje de diversidad (especies) capturado en una muestra por los datos de secuencia actuales. En general, estamos buscando el 95% o mejor. También podemos ver el número de secuencia nseqs y el número de OTU únicassolos en cada muestra.
** summary.single **
Si lo miramos,
vemos que tenemos una muy buena cobertura de estas muestras,> 99%
Normalization
There are two easy options for normalization in mothur. There is no “right” answer for normalization, and it is a rather hot topic in the field. You may find a number of other methods outside of mothur and even papers arguing against normalization at all.
In mothur, you can sub-sample down to your lowest number of sequences in a single sample with sub.sample.
sub.sample
mothur > sub.sample(shared=example.final.an.unique_list.shared)Sampling 12597 from each group.
0.03
Output File Names:
example.final.an.unique_list.0.03.subsample.sharedNote: The default is to subsample to the number of sequences which is seen in the smallest sample. If we has a sample which did not sequence well and only had 10 reads, we would cut out all but 10 reads from every one of our samples using this method. In order to subsample to a specific number, 10000 for example, include “size=10000”.
Or you can calculate percent relative abundance (RA) (# OTU1 sequences/total # sequences in sample) and normalize to your lowest number of sequences in a single sample i.e. %RA * lowest # of sequences
normalize.shared
mothur > normalize.shared(shared=example.final.an.unique_list.shared, method=totalgroup)Normalizing to 12597.
0.03
Output File Names:
example.final.an.unique_list.0.03.norm.sharedNote: The “method” parameter here is used if we simply want to normalize to the number of sequences which is seen in the smallest sample. If we has a sample which did not sequence well and only had 10 reads, we would cut out all but 10 reads from every one of our samples using this method. In order to normalize to a specific number,10000 for example, use “norm=10000” in place of “method=totalgroup”.
When we compare these two methods, we see some differences that can be easily viewed with summary.single
summary.single
mothur > summary.single(shared=example.final.an.unique_list.0.03.subsample.shared, label=0.03, calc=coverage-nseqs)
mothur > summary.single(shared=example.final.an.unique_list.0.03.norm.shared, label=0.03, calc=coverage-nseqs)mothur > system(head -7 example.final.an.unique_list.0.03.subsample.groups.summary)label group coverage nseqs
0.03 L62 0.993887 12597.000000
0.03 S62 0.974121 12597.000000mothur > system(head -7 example.final.an.unique_list.0.03.norm.groups.summary)label group coverage nseqs
0.03 L62 0.993887 12597.000000
0.03 S62 0.971580 12491.000000We see that:
- Sub-sampling takes all samples down to the exact same number of sequences while normalizing takes them near the same number. This occurs in normalization due to rounding, because
.sharedfiles only allow whole numbers. When normalizing to a high number, like >12000 here, this is not an issue as the difference is only ~1%. However, if normalizing to a lower number, this can cause the same problems with unequal sample coverage as we were trying to avoid by normalizing in the first place! In this case, I suggest normalizing outside of mothur in a program, like R, that would allow non-whole numbers in your OTU table. - Our coverage went down for the
S62sample. This is still above our goal of 95%, so no issues. However, if you were to normalize or sub-sample to too low of a number, you will have very poor coverage for more diverse samples. In this case, you should consider removing samples with the lowest number of sequences and normalizing/sub-sampling to a higher value.
Of note, the norm.shared table is the data type used for all R-based analyses in the R workshop.
Normalización Hay dos opciones fáciles para la normalización en mothur. No hay una respuesta “correcta” para la normalización, y es un tema candente en el campo. Puede encontrar otros métodos fuera de mothur e incluso documentos que argumentan en contra de la normalización.
En mothur, puede submuestrear hasta el menor número de secuencias en una sola muestra con sub.sample.
** sub.sample **
** Nota: ** El valor predeterminado es submuestra al número de secuencias que se ve en la muestra más pequeña. Si tenemos una muestra que no se secuenció bien y solo tenía 10 lecturas, recortaríamos todas las lecturas menos 10 de cada una de nuestras muestras con este método. Para submuestrear a un número específico, 10000 por ejemplo, incluya “tamaño = 10000”.
O puede calcular el porcentaje de abundancia relativa (RA) (secuencias # OTU1 / # total de secuencias en la muestra) y normalizar al menor número de secuencias en una sola muestra , es decir, % RA * el menor número de secuencias
** normalize.shared **
** Nota: ** El parámetro “método” aquí se usa si simplemente queremos normalizar el número de secuencias que se ve en la muestra más pequeña. Si tenemos una muestra que no se secuenció bien y solo tenía 10 lecturas, recortaríamos todas las lecturas menos 10 de cada una de nuestras muestras con este método. Para normalizar a un número específico, 10000 por ejemplo, use “norma = 10000” en lugar de “método = grupo total”.
Cuando comparamos estos dos métodos, vemos algunas diferencias que se pueden ver fácilmente con summary.single
** summary.single **
Vemos eso:
- El submuestreo toma todas las muestras hasta el * mismo exacto * número de secuencias mientras que la normalización las toma * cerca * del mismo número. Esto ocurre en la normalización debido al redondeo, porque los archivos
.sharedsolo permiten números enteros. Cuando se normaliza a un número alto, como> 12000 aquí, esto no es un problema ya que la diferencia es solo ~ 1%. Sin embargo, si se normaliza a un número inferior, ¡esto puede causar los mismos problemas con una cobertura de muestra desigual ya que estábamos intentando evitar la normalización en primer lugar! En este caso, sugiero que se normalice fuera de mothur en un programa, como R, que permitiría números no enteros en su tabla de OTU. - Nuestra cobertura bajó para la muestra
S62. Esto sigue por encima de nuestro objetivo del 95%, por lo que no hay problemas. Sin embargo, si tuviera que normalizar o submuestrear a un número demasiado bajo, tendrá una cobertura muy baja para muestras más diversas. En este caso, debería considerar eliminar muestras con el menor número de secuencias y normalizar / submuestrear a un valor más alto.
** De nota, la tabla norm.shared es el tipo de datos utilizado para todos los análisis basados en R en el taller R. **
Basic analysis
There are a number of functions in mothur to perform statistical and visual analysis of your microbiota data. However, mothur itself, cannot display figures and does not handle much manipulation of the data. So, most analyses should be done in a another program, like R, and this section is just to give you some tools to begin to explore your data in mothur.
We will focus on visuals here as statistical modeling is outside the scope of this tutorial and specific to a given data set. Please note that mothur’s figures are not the prettiest…
If you would like to explore statistical analyses in mothur, please see the commands page for:
- Amova: Analysis of molecular variance
- Anosim: Analysis of similarity
- Cooccurence
- Homova: Homogeneity of molecular variance
- Kruskal.wallis
- Otu.association
- Parsimony
- Unifrac, weighted and unweighted
Análisis básico Hay una serie de funciones en mothur para realizar análisis estadísticos y visuales de los datos de su microbiota. Sin embargo, mothur en sí no puede mostrar las cifras y no maneja mucha manipulación de los datos. Entonces, la mayoría de los análisis deben hacerse en otro programa, como R, y esta sección es solo para darle algunas herramientas para comenzar a explorar sus datos en mothur.
Nos enfocaremos en las imágenes aquí ya que el modelado estadístico está fuera del alcance de este tutorial y es específico para un conjunto de datos dado. Tenga en cuenta que las figuras de mothur no son las más bonitas …
Si desea explorar análisis estadísticos en mothur, consulte la [página de comandos] (https://www.mothur.org/wiki/Category:Commands) para:
- Amova: análisis de varianza molecular
- Anosim: análisis de similitud
- Cooccurence
- Homova: Homogeneidad de la varianza molecular
- Kruskal.wallis
- Otu.association
- Parsimonia
- Unifrac, ponderado y no ponderado
Alpha diversity
A lot of research is interested not in the specific microorganisms in a community but the overall diversity and richness of that community. This is within sample diversity. You can calculate all of the classic metrics in mothur. For a full list, see here.
The most important thing before calculating diversity is to normalize your data. No matter the actual community, a sample with 5000 sequences will appear more diverse than one with 3000 due to sequencing errors. This can be done using the norm or subsample data we made above or with the subsample=t modifier, which randomly subsamples a number of times and averages across those subsamples.
We will calculate alpha diversity and richness using summary.single.
summary.single
mothur > summary.single(shared=example.final.an.unique_list.0.03.subsample.shared, label=0.03, calc=coverage-nseqs-bergerparker-chao-shannon)
mothur > summary.single(shared=example.final.an.unique_list.0.03.norm.shared, label=0.03, calc=coverage-nseqs-bergerparker-chao-shannon)Output File Names:
example.final.an.unique_list.0.03.subsample.groups.summary
Output File Names:
example.final.an.unique_list.0.03.norm.groups.summaryOr using subsample=t
mothur > summary.single(shared=example.final.an.unique_list.shared, label=0.03, subsample=t, calc=coverage-nseqs-bergerparker-chao-shannon)Output File Names:
example.final.an.unique_list.groups.ave-std.summary
example.final.an.unique_list.groups.summaryWhen we compare alpha-diversity with the two normalization methods…
mothur > system(head -7 example.final.an.unique_list.0.03.subsample.groups.summary)label group coverage nseqs bergerparker chao chao_lci chao_hci shannon shannon_lci shannon_hci
0.03 L62 0.993887 12597.000000 0.171152 822.357664 812.442451 840.864694 5.151479 5.115362 5.187596
0.03 S62 0.974121 12597.000000 0.173533 1324.876812 1228.290779 1453.934689 4.858460 4.820207 4.896713mothur > system(head -7 example.final.an.unique_list.0.03.norm.groups.summary)label group coverage nseqs bergerparker chao chao_lci chao_hci shannon shannon_lci shannon_hci
0.03 L62 0.993887 12597.000000 0.171152 822.357664 812.442451 840.864694 5.151479 5.115362 5.187596
0.03 S62 0.971580 12491.000000 0.182211 1387.900000 1286.258494 1522.104844 4.826471 4.787373 4.865569mothur > system(head -7 example.final.an.unique_list.groups.summary)label group coverage nseqs bergerparker chao chao_lci chao_hci shannon shannon_lci shannon_hci
0.03 L62 0.993887 12597.000000 0.171152 822.357664 812.442451 840.864694 5.151479 5.115362 5.187596
0.03 S62 0.997046 107972.000000 0.180658 1958.632353 1896.839661 2040.861481 4.898438 4.884562 4.912314we see that richness (chao) is more effected by normalization method than diversity (shannon).
Also, mothur gives us the low confidence interval (lci) and high confidence interval (hci) for our alpha-metrics. These could be used as error bars, if desired.
Alpha diversity Mucha investigación no está interesada en los microorganismos específicos de una comunidad, sino en la diversidad y riqueza general de esa comunidad. Esto es * dentro de la diversidad de la muestra *. Puede calcular todas las métricas clásicas en mothur. Para obtener una lista completa, consulte [aquí] (https://www.mothur.org/wiki/Calculators).
Lo más importante antes de calcular la diversidad es ** normalizar sus datos **. No importa la comunidad real, una muestra con 5000 secuencias aparecerá más diversa que una con 3000 debido a errores de secuencia. Esto se puede hacer usando la norma o los datos de la submuestra que hicimos anteriormente o con el modificador subsample = t, que muestra aleatoriamente un número de veces y promedios a través de esas submuestras.
Calcularemos la diversidad alfa y la riqueza usando summary.single.
** summary.single **
O usando submuestra = t
Cuando comparamos la diversidad alfa con los dos métodos de normalización …
vemos que la riqueza (chao) se ve más afectada por el método de normalización que por la diversidad (shannon).
Además, mothur nos da el intervalo de confianza bajo (lci) y el intervalo de confianza alto (hci) para nuestras métricas alfa. Estos podrían usarse como barras de error, si así lo desea.
Beta diversity
This is between sample diversity meaning every pairwise comparison of two samples has a unique value. You can visualize beta-diversity in a number of ways.
PCA
For principle component analysis (PCA), an overall microbial community is represented by a single point in an xy- or xyz-plane. Two points that are closer to each other indicate that the overall microbiota of those two samples are more similar than two points that are farther apart.
PCA fits axes to the data until 100% of the variation is explained. Thus, you do not specify the number of axes to calculate; you just plot xy or xyz after reviewing how much adding the z-axis adds to the representation of the data. PCA is calculated from the shared file.
pca
mothur > pca(shared=example.final.an.unique_list.0.03.norm.shared, label=0.03)Processing 0.03
Rsq 1 axis: 0
Rsq 2 axis: 0
Rsq 3 axis: 0
Output File Names:
example.final.an.unique_list.0.03.norm.0.03.pca.axes
example.final.an.unique_list.0.03.norm.0.03.pca.loadingsThe outputs are
- axes: the coordinate values for each point (sample) up to N axes
- loadings: the percent of variation in the data explained by each axis
Making the plot (Excel), we have only 2 samples in this tutorial so we only need 1 axis to explain and separate our data. This will not occur with full data sets.
Beta diversity Esta es * entre la diversidad de la muestra * lo que significa que cada comparación por pares de dos muestras tiene un valor único. Puede visualizar la diversidad beta de varias maneras.
PCA Para el análisis de componentes principales (PCA), una comunidad microbiana global está representada por un único punto en un plano xy o xyz. Dos puntos que están más cerca uno del otro indican que la microbiota general de esas dos muestras es más similar que dos puntos que están más separados.
PCA ajusta los ejes a los datos hasta que se explique el 100% de la variación. Por lo tanto, no se especifica el número de ejes para calcular; simplemente traza xy o xyz después de revisar cuánto agrega el eje z a la representación de los datos. PCA se calcula a partir del archivo compartido.
** pca **
Los resultados son
- ejes: los valores de coordenadas para cada punto (muestra) hasta N ejes
- cargas: el porcentaje de variación en los datos explicados por cada eje
Al hacer la trama (Excel), tenemos solo 2 muestras en este tutorial, por lo que solo necesitamos 1 eje para explicar y separar nuestros datos. Esto no ocurrirá con los conjuntos de datos completos.
! [] (PCA.png)
nMDS
nMDS is similar to PCA in that it assigns a single xy(z) point to each sample relative to all other samples. nMDS is more robust than PCA, however, because it permutes through calculations a number of times to find a best fit. This fit is constrained within the number of axes you specify.
mothur currently does not have a function for calculating nMDS for samples. The current nmds function calculates a point for each OTU in a plane. However, this method is something to consider when analyzing data outside of mothur.
nMDS nMDS es similar a PCA en que asigna un único punto xy (z) a cada muestra en relación con todas las demás muestras. nMDS es más robusto que PCA, sin embargo, porque permuta a través de cálculos varias veces para encontrar el mejor ajuste. Este ajuste está restringido dentro del número de ejes que especifique.
mothur actualmente no tiene una función para calcular nMDS para muestras. La función nmds actual calcula un punto para cada OTU en un plano. Sin embargo, este método es algo a considerar cuando se analizan datos fuera de mothur.
Heatmaps of samples
You can visualize beta-diversity metrics by heatmaps with heatmap.sim. Each box is a pairwise comparison between two samples. mothur has several beta-diversity calculator options. Here, we will use the Bray-Curtis metric.
heatmap.sim
mothur > heatmap.sim(shared=example.final.an.unique_list.0.03.norm.shared, label=0.03, calc=braycurtis)0.03
Output File Names:
example.final.an.unique_list.0.03.norm.0.03.braycurtis.heatmap.sim.svg
Mapas de calor de muestras Puede visualizar métricas de diversidad beta mediante heatmaps con heatmap.sim. Cada caja es una comparación por pares entre dos muestras. mothur tiene varias [opciones] de calculadora de diversidad beta (https://www.mothur.org/wiki/Calculators). Aquí, usaremos la métrica de Bray-Curtis.
** heatmap.sim **
Trees of samples
You can easily make a tree where each node is an individual sample in mothur using tree.shared. We will run it with our normalized shared file and Bray-Curtis here.
tree.shared
mothur > tree.shared(shared=example.final.an.unique_list.0.03.norm.shared, label=0.03, calc=braycurtis)Using 2 processors.
0.03
Output File Names:
example.final.an.unique_list.0.03.norm.braycurtis.0.03.treIf we open it in a program such as FigTree: 
Árboles de muestras Puede hacer fácilmente un árbol donde cada nodo sea una muestra individual en mothur usando tree.shared. Lo ejecutaremos con nuestro archivo “compartido” normalizado y Bray-Curtis aquí.
** tree.shared **
Si lo abrimos en un programa como [FigTree] (http://tree.bio.ed.ac.uk/software/figtree/): ! [] (tree_shared.png)
OTU analyses
Trees of OTUs
mothur does not have an easy way to create a tree where the nodes are OTUs. In order to make a tree of your OTUs outside of mothur, you need a representative sequence for each OTU. Otherwise, you end up with many more branches then there are OTUs (i.e. species) in your data set. You then calculate a distance matrix of these representative OTU sequences and calculate a tree from this in R or another program. This sort of tree is used for UniFrac beta-diversity calculators.
We will not actually run this analysis in this workshop as it is RAM and time intensive. However, if we were to extract representative sequences, we would use get.oturep. The first option is to have mothur calculate the best possible representative sequence for each OTU. This is the sequence with smallest maximum distance to the other sequences in an OTU group.
DO NOT RUN THIS
mothur > get.oturep(column=example.final.dist, list=example.final.an.unique_list.list, fasta=example.final.fasta, count=example.final.count_table, label=0.03, method=distance, large=T)Or you can have mothur just choose the most abundant unique sequence within each OTU group as the representative. This representative sequence may not be the best possible one for all uniques within an OTU, but this runs much faster and with much less RAM than the other method.
DO NOT RUN THIS
mothur > get.oturep(list=example.final.an.unique_list.list, fasta=example.final.fasta, count=example.final.count_table, label=0.03, method=abundance)These representative sequences have long, horrible names given to them by the MiSeq. We rename them based on the OTU# using a custom python script, courtesy of Tony Neumann (Suen lab).
import sys
import re
#an empty list to hold the lines, as strings, of the input fasta file
info = []
#iterate through the input file and file up the list
with open(sys.argv[1], 'r') as fasta_in:
for line in fasta_in:
info.append(line.replace('\r', '').rstrip('\n'))
#an empty list to hold the OTU fasta seqs and their IDs as tuples
seqs = []
#iterate through the info list, replace the complex seq IDs with just the OTU ID, and fill up the seqs list
for index, item in enumerate(info):
if index == 0:
ID = re.search('Otu[0-9]*',item[1:]).group(0)
else:
if index % 2 == 0:
ID = re.search('Otu[0-9]*',item[1:]).group(0)
else:
seqs.append((ID, item))
#write the contents of the seqs list to a new file called "clean_repFasta.fasta"
with open("clean_repFasta.fasta", 'w') as fasta_out:
for item in seqs:
fasta_out.write('>' + item[0] + '\n' + item[1] + '\n')This outputs a new fasta called clean_repFasta.fasta. Then, we would calculate the distance matrix with the new fasta. We definitely what to use the lower triangle lt output here so our matrix is small enough to read into R later on.
DO NOT RUN THIS
mothur > dist.seqs(fasta=clean_repFasta.fasta, output=lt)You would then load this distance matrix into R or another program to calculate and visualize the tree.
OTU analiza
Árboles de OTUs mothur no tiene una manera fácil de crear un árbol donde los nodos son OTU. Para hacer un árbol de tus OTU fuera de mothur, necesitas una secuencia representativa para cada OTU. De lo contrario, terminas con muchas más ramas, entonces hay OTU (* i.e. * species) en tu conjunto de datos. A continuación, calcula una matriz de distancia de estas secuencias OTU representativas y calcula un árbol a partir de esto en R u otro programa. Este tipo de árbol se usa para las calculadoras * UniFrac beta-diversity *.
En realidad, no ejecutaremos este análisis en este taller, ya que es RAM y requiere mucho tiempo. Sin embargo, si tuviéramos que extraer secuencias representativas, usaríamos get.oturep. La primera opción es hacer que mothur calcule la mejor secuencia representativa posible para cada OTU. Esta es la secuencia con la distancia máxima más pequeña a las otras secuencias en un grupo OTU.
** NO EJECUTE ESTO **
O puede tener mothur simplemente elija la secuencia única más abundante dentro de cada grupo OTU como representante. Esta secuencia representativa puede no ser la mejor posible para todos los únicos dentro de una OTU, pero esto se ejecuta mucho más rápido y con mucha menos memoria RAM que el otro método.
** NO EJECUTE ESTO **
Estas secuencias representativas tienen nombres largos y horribles que les da el MiSeq. Los cambiamos de nombre en función de la OTU # utilizando una secuencia de comandos personalizada de python, cortesía de Tony Neumann (laboratorio de Suen).
Esto produce un nuevo fasta llamado clean_repFasta.fasta. Entonces, calcularíamos la matriz de distancia con el nuevo fasta. Definitivamente usamos el resultado del triángulo inferior lt aquí, así que nuestra matriz es lo suficientemente pequeña para leer en R más adelante.
** NO EJECUTE ESTO **
Luego cargaría esta matriz de distancia en R u otro programa para calcular y visualizar el árbol.
Heatmaps of OTUs
Similar to the heatmap of samples, we can do a heatmap for OTUs. Instead of pairwise comparisons, each box represents the relative abundance of an OTU in a sample. You can scale the abundance to make differences and low abundance OTUs easier to see. We will use the default log10 scale and only visualize the 10 most abundant OTUs.
heatmap.bin
mothur > heatmap.bin(shared=example.final.an.unique_list.0.03.norm.shared, label=0.03, scale=log10, numotu=10)0.03
Output File Names:
example.final.an.unique_list.0.03.norm.0.03.heatmap.bin.svg
Mapas de calor de OTUs Similar al mapa de calor de las muestras, podemos hacer un mapa de calor para OTU. En lugar de comparaciones por pares, cada cuadro representa la abundancia relativa de una OTU en una muestra. Puede escalar la abundancia para hacer que las diferencias y la baja abundancia de OTU sean más fáciles de ver. Utilizaremos la escala log10 predeterminada y solo visualizaremos las 10 OTU más abundantes.
** heatmap.bin **
Venns
Venn diagrams can show you which OTUs are shared by samples. mothur can do 2, 3, or 4 sample Venns.
Note: If you have more than 4 samples, you can use the permute parameter to have mothur create all Venn groups of 2, 3, or 4 between all of your samples.
venn
mothur > venn(shared=example.final.an.unique_list.0.03.norm.shared, label=0.03, calc=sharedsobs)0.03
Output File Names:
example.final.an.unique_list.0.03.norm.0.03.sharedsobs.L62-S62.svg
example.final.an.unique_list.0.03.norm.0.03.sharedsobs.L62-S62.sharedotus
Venns Los diagramas de Venn pueden mostrarle qué OTU comparten las muestras. mothur puede hacer 2, 3 o 4 muestras de Venns.
** Nota **: si tiene más de 4 muestras, puede usar el parámetro permute para que mothur cree todos los grupos de Venn de 2, 3 o 4 entre todas sus muestras.
** venn **
So your distance matrix is too large…
Large data sets (many samples) and those with high error (V3-4) can create very large .dist files on the order of terabytes. This can prevent downstream analysis like clustering because the distance matrix is too large to read into your available RAM. And even when things do run, you may not want to wait days for a step to complete. If this happens, there are a number of steps you can take to reduce the distance matrix.
Solutions that do not remove sequences
dist.seqs
If you’re only going to look at genus (0.05) and/or species (0.03) OTUs in the end, set cutoff=0.05 for dist.seqs so it doesn’t save values larger than that.
You can also save the output in lower triangle format so that only one value is saved for every pairwise comparison between sequences. This is done with output=lt in dist.seqs. If you do this, the output will be example.final.phylip.dist and you will need to alter cluster.split to be phylip=example.final.phylip.dist instead of using column=.
cluster.split
Run cluster with large=true to tell mothur the distance matrix is too large to read into RAM all at once.
pre.cluster
If you have a longer amplicon, less overlap of the two sequences that are combined into contigs, and/or a high error rate in sequencing, it is worth increasing pre.cluster to diffs=3 or 4. This will group more sequences together and while it does not remove any data, you will end up with fewer uniques and therefore, a smaller distance matrix.
Solutions that do remove sequences
split.abund
This function can be used to remove very low abundance sequences. It creates a rare and abund data set. When cutoff=1, the abund data set includes all sequences that are present in the overall data set in 2 or more copies. The rare data set is often referred to as singletons. Removing singletons (or doubletons, etc) dramatically reduces your unique sequences without strongly impacting beta-diversity. However, alpha-diversity is effected by singletons so please be careful to compare values only across data sets that have been handled in a similar way. This command would be used as:
mothur > split.abund(fasta=example.final.fasta, count=example.final.count_table, cutoff=1)
Entonces tu matriz de distancia es demasiado grande … Los grandes conjuntos de datos (muchas muestras) y aquellos con alto error (V3-4) pueden crear archivos .dist muy grandes del orden de los terabytes. Esto puede evitar el análisis posterior como la agrupación porque la matriz de distancia es demasiado grande para leer en la memoria RAM disponible. E incluso cuando las cosas se ejecutan, es posible que no desee esperar días para completar un paso. Si esto sucede, hay una serie de pasos que puede seguir para reducir la matriz de distancia.
Soluciones que ** no ** eliminan secuencias
** dist.seqs **
Si solo vas a mirar OTUs genus (0.05) y / o species (0.03) al final, establece cutoff = 0.05 para dist.seqs para que no guarde valores más grandes que eso.
También puede guardar la salida en formato de triángulo inferior para que solo se guarde un valor por cada comparación por pares entre secuencias. Esto se hace con output = lt en dist.seqs. Si haces esto, la salida será example.final.phylip.dist y necesitarás alterar cluster.split para que seaphylip = example.final.phylip.dist en lugar de usar column =.
** cluster.split **
Ejecute el clúster con large = true para decirle a mothur que la matriz de distancia es demasiado grande para leer en la RAM de una sola vez.
** pre.cluster **
Si tienes un amplicón más largo, menos superposición de las dos secuencias que se combinan en contigs, y / o una alta tasa de error en la secuencia, vale la pena aumentar pre.cluster a diffs = 3 o 4. Esto agrupará más secuencias juntas y aunque no elimina ningún dato, terminará con menos únicos y, por lo tanto, una matriz de distancia más pequeña.
Soluciones que ** eliminan ** secuencias
** split.abund **
Esta función se puede usar para eliminar secuencias de muy baja abundancia. Crea un conjunto de datos raro yabund. Cuando cutoff = 1, el conjunto de datosabund incluye todas las secuencias que están presentes en el conjunto de datos global en 2 o más copias. El conjunto de datos ‘raro’ a menudo se conoce como singletons. La eliminación de singletons (o doubletons, etc.) reduce drásticamente sus secuencias únicas sin afectar fuertemente la diversidad beta. Sin embargo, la diversidad alfa se efectúa por singleton, así que tenga cuidado de comparar los valores solo en los conjuntos de datos que se han manejado de forma similar. Este comando se usaría como:
Batch mode
Now that you know what you are doing, and why, how do you avoid ever having to type all of these ever again? mothur is able to accept batch files which contain a list of all of the commands you want to feed it. Here is a link to the Schloss lab SOP. The zip file with their data includes a batch file which you could adapt: https://www.mothur.org/wiki/MiSeq_SOP
First, let’s get out of mothur…
mothur > quit()… then, let’s copy over our mock batch file.
$ cp /home/capacitaciones/bioinformatica-UCR/batch.txt ~/MiSeq_SOPNow, let’s look at the batch file. It contains three commands which should run fairly quickly: summary.seqs, screen.seqs, and unique.seqs.
$ head -5 batch.txtsummary.seqs(fasta=example.trim.contigs.fasta)
screen.seqs(fasta=example.trim.contigs.fasta, group=example.contigs.groups, maxambig=0, maxlength=300, maxhomop=8)
unique.seqs(fasta=example.trim.contigs.good.fasta)We feed it to mothur using the filepath to the mothur executable followed by a space, then the name of our batch file.
$ module load mothur/1.39.5
$ /opt/bioinf/mothur-1.39.5/mothur batch.txtYou can see that mothur was opened, the commands were executed, then mothur was automatically closed. This is useful for large data sets, where certain commands (pre.cluster, dist.seqs) can take a VERY long time. This allows you to execute all the necessary commands without having to babysit the computer.
Note: Because we are in a new session of mothur, it has forgotten our preferences for multiple processors. You must indicate at the beginning of each session how many processors you want mothur to use. This is particularly important when you are using a batch file with a large data set and not paying as close attention to your data as it runs. Your first command in a batch file should indicate the number of processors you want to use.
Por lotes Ahora que sabes lo que estás haciendo, y por qué, ¿cómo evitarás tener que volver a escribir todo esto alguna vez? mothur puede aceptar archivos por lotes que contienen una lista de todos los comandos que desea alimentar. Aquí hay un enlace al SOP del laboratorio Schloss. El archivo comprimido con sus datos incluye un archivo por lotes que puede adaptar: https://www.mothur.org/wiki/MiSeq_SOP
Primero, salgamos de mothur …
… entonces, copiemos nuestro archivo de lote falso.
Ahora, veamos el archivo por lotes. Contiene tres comandos que deben ejecutarse con bastante rapidez: summary.seqs, screen.seqs y unique.seqs.
Lo alimentamos a mothur utilizando el camino de archivo al ejecutable mothur seguido de un espacio, luego el nombre de nuestro archivo por lotes.
Se puede ver que mothur se abrió, se ejecutaron los comandos, luego mothur se cerró automáticamente. Esto es útil para grandes conjuntos de datos, donde ciertos comandos (pre.cluster, dist.seqs) pueden demorar MUCHO tiempo. Esto le permite ejecutar todos los comandos necesarios sin tener que cuidar a la computadora.
** Nota: ** Debido a que estamos en una nueva sesión de mothur, ha olvidado nuestras preferencias para múltiples procesadores. Debe indicar al comienzo de ** cada sesión ** cuántos procesadores quiere que use Mothur. Esto es particularmente importante cuando utiliza un archivo por lotes con un gran conjunto de datos y no presta tanta atención a sus datos mientras se ejecuta. Su primer comando en un archivo por lotes debe indicar la cantidad de procesadores que desea usar.
Ending session
It is important that we close out of the session at the end. Do so by typing “exit” in the terminal.
$ exit