diff --git a/episodes/04-redirection.md b/episodes/04-redirection.md index a8596b46..5580d877 100644 --- a/episodes/04-redirection.md +++ b/episodes/04-redirection.md @@ -31,134 +31,66 @@ regular expressions in this lesson, and are instead going to specify the strings we are searching for. Let's give it a try! -::::::::::::::::::::::::::::::::::::::::: callout - -## Nucleotide abbreviations - -The four nucleotides that appear in DNA are abbreviated `A`, `C`, `T` and `G`. -Unknown nucleotides are represented with the letter `N`. An `N` appearing -in a sequencing file represents a position where the sequencing machine was not able to -confidently determine the nucleotide in that position. You can think of an `N` as being aNy -nucleotide at that position in the DNA sequence. - -:::::::::::::::::::::::::::::::::::::::::::::::::: - -We'll search for strings inside of our fastq files. Let's first make sure we are in the correct +We'll search for strings inside of a metadata file. Let's first make sure we are in the correct directory: ```bash -$ cd ~/shell_data/untrimmed_fastq +$ cd ~/shell_data/sra_metadata ``` -Let's look for lines that contain `ACGT`. +Let's look for lines that contain `PAIRED`. ```bash -$ grep ACGT SRR098026.fastq +$ grep PAIRED SraRunTable.txt ``` -To get only the number of lines with `ACGT`, we can use the `-c` flag. +To get only the number of lines with `PAIRED`, we can use the `-c` flag. This is useful if you are unsure about the number of lines that will be found. ```bash -$ grep -c ACGT SRR098026.fastq +$ grep -c PAIRED SraRunTable.txt ``` -Suppose we want to see how many reads in our file have really bad segments containing 10 consecutive unknown nucleotides (Ns). - -::::::::::::::::::::::::::::::::::::::::: callout - -## Determining quality - -In this lesson, we're going to be manually searching for strings of `N`s within our sequence -results to illustrate some principles of file searching. It can be really useful to do this -type of searching to get a feel for the quality of your sequencing results, however, in your -research you will most likely use a bioinformatics tool that has a built-in program for -filtering out low-quality reads. You'll learn how to use one such tool in -[a later lesson](https://datacarpentry.org/wrangling-genomics/02-quality-control). - -:::::::::::::::::::::::::::::::::::::::::::::::::: - -Let's search for the string NNNNNNNNNN in the SRR098026 file: +You can use case-insensitive searching with the `-i` flag. +This is useful if you are unsure if what you are searching for is in upper- or lower-case or a mix. ```bash -$ grep NNNNNNNNNN SRR098026.fastq +$ grep -i paired SraRunTable.txt ``` -This command returns a lot of output to the terminal. Every single line in the SRR098026 -file that contains at least 10 consecutive Ns is printed to the terminal, regardless of how long or short the file is. -We may be interested not only in the actual sequence which contains this string, but -in the name (or identifier) of that sequence. We discussed in a previous lesson -that the identifier line immediately precedes the nucleotide sequence for each read -in a FASTQ file. We may also want to inspect the quality scores associated with -each of these reads. To get all of this information, we will return the line -immediately before each match and the two lines immediately after each match. - -We can use the `-B` argument for grep to return a specific number of lines before -each match. The `-A` argument returns a specific number of lines after each matching line. Here we want the line *before* and the two lines *after* each -matching line, so we add `-B1 -A2` to our grep command: +The `-v` option for `grep` search stands for `--invert-match` meaning `grep` will now only display the +lines which do not match the searched pattern. ```bash -$ grep -B1 -A2 NNNNNNNNNN SRR098026.fastq +$ grep -v SINGLE SraRunTable.txt ``` -One of the sets of lines returned by this command is: +Notice that you now get the header line and the paired-end samples, because these do not match the +pattern `SINGLE`. -```output -@SRR098026.177 HWUSI-EAS1599_1:2:1:1:2025 length=35 -CNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN -+SRR098026.177 HWUSI-EAS1599_1:2:1:1:2025 length=35 -#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -``` ::::::::::::::::::::::::::::::::::::::: challenge ## Exercise -1. Search for the sequence `GNATNACCACTTCC` in the `SRR098026.fastq` file. - Have your search return all matching lines and the name (or identifier) for each sequence - that contains a match. +1. Count the number of single-end samples using a case-insensitive search. -2. Search for the sequence `AAGTT` in both FASTQ files. - Have your search return all matching lines and the name (or identifier) for each sequence - that contains a match. +2. Count the number of single-end samples and the header line. ::::::::::::::: solution ## Solution -1. `grep -B1 GNATNACCACTTCC SRR098026.fastq` +1. `grep -c -i Single SraRunTable.txt` ``` -@SRR098026.245 HWUSI-EAS1599_1:2:1:2:801 length=35 -GNATNACCACTTCCAGTGCTGANNNNNNNGGGATG +35 ``` -2. `grep -B1 AAGTT *.fastq` +2. `grep -c -v PAIRED SraRunTable.txt` ``` -SRR097977.fastq-@SRR097977.11 209DTAAXX_Lenski2_1_7:8:3:247:351 length=36 -SRR097977.fastq:GATTGCTTTAATGAAAAAGTCATATAAGTTGCCATG --- -SRR097977.fastq-@SRR097977.67 209DTAAXX_Lenski2_1_7:8:3:544:566 length=36 -SRR097977.fastq:TTGTCCACGCTTTTCTATGTAAAGTTTATTTGCTTT --- -SRR097977.fastq-@SRR097977.68 209DTAAXX_Lenski2_1_7:8:3:724:110 length=36 -SRR097977.fastq:TGAAGCCTGCTTTTTTATACTAAGTTTGCATTATAA --- -SRR097977.fastq-@SRR097977.80 209DTAAXX_Lenski2_1_7:8:3:258:281 length=36 -SRR097977.fastq:GTGGCGCTGCTGCATAAGTTGGGTTATCAGGTCGTT --- -SRR097977.fastq-@SRR097977.92 209DTAAXX_Lenski2_1_7:8:3:353:318 length=36 -SRR097977.fastq:GGCAAAATGGTCCTCCAGCCAGGCCAGAAGCAAGTT --- -SRR097977.fastq-@SRR097977.139 209DTAAXX_Lenski2_1_7:8:3:703:655 length=36 -SRR097977.fastq:TTTATTTGTAAAGTTTTGTTGAAATAAGGGTTGTAA --- -SRR097977.fastq-@SRR097977.238 209DTAAXX_Lenski2_1_7:8:3:592:919 length=36 -SRR097977.fastq:TTCTTACCATCCTGAAGTTTTTTCATCTTCCCTGAT --- -SRR098026.fastq-@SRR098026.158 HWUSI-EAS1599_1:2:1:1:1505 length=35 -SRR098026.fastq:GNNNNNNNNCAAAGTTGATCNNNNNNNNNTGTGCG +36 ``` ::::::::::::::::::::::::: @@ -179,91 +111,53 @@ use other commands to analyze this data. The command for redirecting output to a file is `>`. -Let's try out this command and copy all the records (including all four lines of each record) -in our FASTQ files that contain -'NNNNNNNNNN' to another file called `bad_reads.txt`. +Let's search for the metadata for sample `SRR097977` and redirect the output to a file. ```bash -$ grep -B1 -A2 NNNNNNNNNN SRR098026.fastq > bad_reads.txt +$ grep SRR097977 SraRunTable.txt > metadata.txt ``` -::::::::::::::::::::::::::::::::::::::::: callout - -## File extensions - -You might be confused about why we're naming our output file with a `.txt` extension. After all, -it will be holding FASTQ formatted data that we're extracting from our FASTQ files. Won't it -also be a FASTQ file? The answer is, yes - it will be a FASTQ file and it would make sense to -name it with a `.fastq` extension. However, using a `.fastq` extension will lead us to problems -when we move to using wildcards later in this episode. We'll point out where this becomes -important. For now, it's good that you're thinking about file extensions! - -:::::::::::::::::::::::::::::::::::::::::::::::::: - -The prompt should sit there a little bit, and then it should look like nothing -happened. But type `ls`. You should see a new file called `bad_reads.txt`. +Type `ls`. You should see a new file called `metadata.txt`. We can check the number of lines in our new file using a command called `wc`. `wc` stands for **word count**. This command counts the number of words, lines, and characters -in a file. The FASTQ file may change over time, so given the potential for updates, -make sure your file matches your instructor's output. - -As of Sept. 2020, wc gives the following output: +in a file. ```bash -$ wc bad_reads.txt +$ wc metadata.txt ``` ```output - 802 1338 24012 bad_reads.txt + 1 31 228 metadata.txt ``` This will tell us the number of lines, words and characters in the file. If we want only the number of lines, we can use the `-l` flag for `lines`. ```bash -$ wc -l bad_reads.txt +$ wc -l metadata.txt ``` ```output -802 bad_reads.txt +1 metadata.txt ``` ::::::::::::::::::::::::::::::::::::::: challenge ## Exercise -How many sequences are there in `SRR098026.fastq`? Remember that every sequence is formed by four lines. +How many entries are there in `SraRunTable.txt`? ::::::::::::::: solution ## Solution ```bash -$ wc -l SRR098026.fastq +$ wc -l SraRunTable.txt ``` ```output -996 -``` - -Now you can divide this number by four to get the number of sequences in your fastq file. - -This can be done using [shell integer arithmetic](https://www.gnu.org/software/bash/manual/html_node/Shell-Arithmetic.html) - -```bash -$ echo $((996/4)) -``` - -Note, this will do integer division - if you need floating point arithmetic you can use [bc - an arbitrary precision calculator](https://www.gnu.org/software/bc/manual/html_mono/bc.html) - - -```bash -$ echo "996/4" | bc -``` - -```output -249 +37 SraRunTable.txt ``` ::::::::::::::::::::::::: @@ -274,116 +168,63 @@ $ echo "996/4" | bc ## Exercise -How many sequences in `SRR098026.fastq` contain at least 3 consecutive Ns? +How many paired-end read samples are there in `SraRunTable.txt`? These samples will have metadata that contains the keyword `PAIRED`. ::::::::::::::: solution ## Solution ```bash -$ grep NNN SRR098026.fastq > bad_reads.txt -$ wc -l bad_reads.txt +$ grep PAIRED SraRunTable.txt > metadata.txt +$ wc -l metadata.txt ``` ```output -249 +2 metadata.txt ``` ::::::::::::::::::::::::: :::::::::::::::::::::::::::::::::::::::::::::::::: -We might want to search multiple FASTQ files for sequences that match our search pattern. +We might want to search our file for multiple patterns, e.g. all single-end and all paired-end samples. However, we need to be careful, because each time we use the `>` command to redirect output to a file, the new output will replace the output that was already present in the file. This is called "overwriting" and, just like you don't want to overwrite your video recording of your kid's first birthday party, you also want to avoid overwriting your data files. -```bash -$ grep -B1 -A2 NNNNNNNNNN SRR098026.fastq > bad_reads.txt -$ wc -l bad_reads.txt -``` - -```output -802 bad_reads.txt -``` - -```bash -$ grep -B1 -A2 NNNNNNNNNN SRR097977.fastq > bad_reads.txt -$ wc -l bad_reads.txt -``` - -```output -0 bad_reads.txt -``` - -Here, the output of our second call to `wc` shows that we no longer have any lines in our `bad_reads.txt` file. This is -because the second file we searched (`SRR097977.fastq`) does not contain any lines that match our -search sequence. So our file was overwritten and is now empty. - -We can avoid overwriting our files by using the command `>>`. `>>` is known as the "append redirect" and will -append new output to the end of a file, rather than overwriting it. +Find the paired-end samples in the `SraRunTable.txt` file and take a look at the output with `less`. +Remember you can exit `less` by pressing `q`. ```bash -$ grep -B1 -A2 NNNNNNNNNN SRR098026.fastq > bad_reads.txt -$ wc -l bad_reads.txt +$ grep PAIRED SraRunTable.txt > metadata.txt +$ less metadata.txt ``` -```output -802 bad_reads.txt -``` +Find the single-end samples in the `SraRunTable.txt` file. ```bash -$ grep -B1 -A2 NNNNNNNNNN SRR097977.fastq >> bad_reads.txt -$ wc -l bad_reads.txt +$ grep SINGLE SraRunTable.txt > metadata.txt +$ less metadata.txt ``` -```output -802 bad_reads.txt -``` - -The output of our second call to `wc` shows that we have not overwritten our original data. - -We can also do this with a single line of code by using a wildcard: - -```bash -$ grep -B1 -A2 NNNNNNNNNN *.fastq > bad_reads.txt -$ wc -l bad_reads.txt -``` - -```output -802 bad_reads.txt -``` - -::::::::::::::::::::::::::::::::::::::::: callout - -## File extensions - part 2 +Notice that the paired-end samples are no longer present in the output `metadata.txt` file. This is +because the our second search overwrote the results of the first search. -This is where we would have trouble if we were naming our output file with a `.fastq` extension. -If we already had a file called `bad_reads.fastq` (from our previous `grep` practice) -and then ran the command above using a `.fastq` extension instead of a `.txt` extension, `grep` -would give us a warning. +We can avoid overwriting our files by using the command `>>`. `>>` is known as the "append redirect" and will append new output to the end of a file, rather than overwriting it. ```bash -grep -B1 -A2 NNNNNNNNNN *.fastq > bad_reads.fastq -``` - -```output -grep: input file ‘bad_reads.fastq' is also the output +$ grep PAIRED SraRunTable.txt > metadata.txt +$ grep SINGLE SraRunTable.txt >> metadata.txt +$ less metadata.txt ``` -`grep` is letting you know that the output file `bad_reads.fastq` is also included in your -`grep` call because it matches the `*.fastq` pattern. Be careful with this as it can lead to -some unintended results. - -:::::::::::::::::::::::::::::::::::::::::::::::::: +Note that the paired-end samples are the first two lines of the file and the single-end samples come after (appended). Since we might have multiple different criteria we want to search for, -creating a new output file each time has the potential to clutter up our workspace. We also -thus far haven't been interested in the actual contents of those files, only in the number of -reads that we've found. We created the files to store the reads and then counted the lines in -the file to see how many reads matched our criteria. There's a way to do this, however, that -doesn't require us to create these intermediate files - the pipe command (`|`). +creating a new output file each time has the potential to clutter up our workspace. +We've been redirecting output to a file and then using `less` to view the contents. +There's a way to do this that doesn't require us to create these intermediate files - the pipe command (`|`). This is probably not a key on your keyboard you use very much, so let's all take a minute to find that key. @@ -397,71 +238,30 @@ look at it, like we can with `less`. Well it turns out that we can! We can redir from our `grep` call through the `less` command. ```bash -$ grep -B1 -A2 NNNNNNNNNN SRR098026.fastq | less +$ grep SINGLE SraRunTable.txt | less ``` We can now see the output from our `grep` call within the `less` interface. We can use the up and down arrows to scroll through the output and use `q` to exit `less`. -If we don't want to create a file before counting lines of output from our `grep` search, we could directly pipe -the output of the grep search to the command `wc -l`. This can be helpful for investigating your output if you are not sure -you would like to save it to a file. +::::::::::::::::::::::::::::::::::::::::: callout -```bash -$ grep -B1 -A2 NNNNNNNNNN SRR098026.fastq | wc -l -``` +## Viewing files that are too wide for the terminal -Because we asked `grep` for all four lines of each FASTQ record, we need to divide the output by -four to get the number of sequences that match our search pattern. Since 802 / 4 = 200.5 and we -are expecting an integer number of records, there is something added or missing in `bad_reads.txt`. -If we explore `bad_reads.txt` using `less`, we might be able to notice what is causing the uneven -number of lines. Luckily, this issue happens by the end of the file so we can also spot it with `tail`. +`less` will wrap lines in your terminal if they are too long to be displayed. +Use `less -S` to avoid line-wrapping. +You can use the left and right arrows to scroll across the output similarly to up and down. -```bash -$ grep -B1 -A2 NNNNNNNNNN SRR098026.fastq > bad_reads.txt -$ tail bad_reads.txt -``` +:::::::::::::::::::::::::::::::::::::::::::::::::: -```output -@SRR098026.133 HWUSI-EAS1599_1:2:1:0:1978 length=35 -ANNNNNNNNNTTCAGCGACTNNNNNNNNNNGTNGN -+SRR098026.133 HWUSI-EAS1599_1:2:1:0:1978 length=35 -#!!!!!!!!!##########!!!!!!!!!!##!#! --- --- -@SRR098026.177 HWUSI-EAS1599_1:2:1:1:2025 length=35 -CNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN -+SRR098026.177 HWUSI-EAS1599_1:2:1:1:2025 length=35 -#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -``` - -The fifth and six lines in the output display "--" which is the default action for `grep` to separate groups of -lines matching the pattern, and indicate groups of lines which did not match the pattern so are not displayed. -To fix this issue, we can redirect the output of grep to a second instance of `grep` as follows. +If we don't want to create a file before counting lines of output from our `grep` search, we could directly pipe +the output of the grep search to the command `wc -l`. This can be helpful for investigating your output if you are not sure +you would like to save it to a file. ```bash -$ grep -B1 -A2 NNNNNNNNNN SRR098026.fastq | grep -v '^--' > bad_reads.fastq -$ tail bad_reads.fastq +$ grep SINGLE SraRunTable.txt | wc -l ``` -```output -+SRR098026.132 HWUSI-EAS1599_1:2:1:0:320 length=35 -#!!!!!!!!!##########!!!!!!!!!!##!#! -@SRR098026.133 HWUSI-EAS1599_1:2:1:0:1978 length=35 -ANNNNNNNNNTTCAGCGACTNNNNNNNNNNGTNGN -+SRR098026.133 HWUSI-EAS1599_1:2:1:0:1978 length=35 -#!!!!!!!!!##########!!!!!!!!!!##!#! -@SRR098026.177 HWUSI-EAS1599_1:2:1:1:2025 length=35 -CNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN -+SRR098026.177 HWUSI-EAS1599_1:2:1:1:2025 length=35 -#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -``` - -The `-v` option in the second `grep` search stands for `--invert-match` meaning `grep` will now only display the -lines which do not match the searched pattern, in this case `'^--'`. The caret (`^`) is an **anchoring** -character matching the beginning of the line, and the pattern has to be enclose by single quotes so `grep` does -not interpret the pattern as an extended option (starting with --). - ::::::::::::::::::::::::::::::::::::::::: callout ## Custom `grep` control @@ -530,7 +330,7 @@ foo is abcEFG Let's write a for loop to show us the first two lines of the fastq files we downloaded earlier. You will notice the shell prompt changes from `$` to `>` and back again as we were typing in our loop. The second prompt, `>`, is different to remind us that we haven't finished typing a complete command yet. A semicolon, `;`, can be used to separate two commands written on a single line. ```bash -$ cd ../untrimmed_fastq/ +$ cd ~/shell_data/untrimmed_fastq/ ``` ```bash diff --git a/episodes/05-writing-scripts.md b/episodes/05-writing-scripts.md index f8bbf66c..0da2045a 100644 --- a/episodes/05-writing-scripts.md +++ b/episodes/05-writing-scripts.md @@ -152,6 +152,8 @@ $ nano bad-reads-script.sh Bad reads have a lot of N's, so we're going to look for `NNNNNNNNNN` with `grep`. We want the whole FASTQ record, so we're also going to get the one line above the sequence and the two lines below. We also want to look in all the files that end with `.fastq`, so we're going to use the `*` wildcard. +`--` is the default action for `grep` to separate groups of lines matching the pattern. The caret (`^`) is an **anchoring** character matching the beginning of the line, and the pattern has to be enclosed by single quotes so `grep` does not interpret the pattern as an extended option (starting with --). + ```bash grep -B1 -A2 -h NNNNNNNNNN *.fastq | grep -v '^--' > scripted_bad_reads.txt ```