--- /dev/null
+<style type="text/css">
+.answer {
+display: none;
+margin-top: 10px;
+margin-bottom: 10px;
+margin-left: 30px;
+margin-right: 30px;
+padding: 10px;ul>
+border: 1px dotted grey;
+background-color:#FFFFF2;
+}
+ul > li {
+ padding: 10px 0;
+}
+ul {
+ margin-left:0;
+ list-style-type: none;
+}
+</style>
+
+<script type="text/javascript">// <![CDATA[
+function toggle(id) {
+ e = document.getElementById(id);
+ if (e.style.display == "block") {
+ e.style.display = "none";
+ } else {
+ e.style.display = "block";
+ }
+}
+// ]]></script>
+
+
+Let's try out your new shell skills on some real data.
+
+The following file is a small sample (1%) of a very large text file containing human genetics data. Specifically, it describes genetic variation in three African individuals sequenced as part of the <a href="http://www.1000genomes.org">1000 Genomes Project</a>. Please download it from the following url, open up a shell, and go into the directory where it was saved:
+
+<a href="http://svn.software-carpentry.org/swc/data/1000gp.vcf">http://svn.software-carpentry.org/swc/data/1000gp.vcf</a>
+
+
+<h2>Exercise Part 1 (setup)</h2>
+
+<ul>
+<li>
+If you had forgotten where you downloaded the file, how would you locate the path of all files with that name on the computer (using the shell)?
+(<a href="javascript:toggle('q1_h');">Click for a hint</a>)
+<div id="q1_h" class="answer">
+<code>$ man find</code>
+</div>
+(<a href="javascript:toggle('q1_a');">Click for our answer</a>)
+<div id="q1_a" class="answer">
+<code>$ find / -name "1000gp.vcf"</code>
+</div>
+</li>
+
+
+<li>
+It's usually a good idea to use an empty directory as a workspace so that other files don't get in the way (or accidentally get overwritten or deleted). Create a new subdirectory directory named "sandbox", move our data file there, and make the directory your current working directory (sandbox should be the last part of the path given when you type <code>pwd</code>.
+
+(<a href="javascript:toggle('q2_a');">Click for our answer</a>)
+<div id="q2_a" class="answer">
+<code>$ mkdir sandbox
+$ mv /home/orion/Downloads/1000gp.vcf sandbox
+$ cd sandbox
+</code>
+</div>
+</li>
+</ul>
+
+
+<h2>Exercise Part 2 (analysis)</h2>
+<ul>
+<li>
+The data file you downloaded is a line-based text file. The "vcf" extension lets us know that it's in a specific text format, namely "Variant Call Format". The file starts with a bunch of comment lines (they start with "#" or "##"), and then a large number of data lines. The human genome can be thought of as an encyclopedia, where each chromosome is a volume. Each volume is just a long string of characters, but rather than the english alphabet, the genome uses just the characters "A", "C", "G", and "T". This VCF file lists the differences between the three African individuals and a standard "individual" called the reference (actually based upon a few different people). Each line in the file corresponds to a difference. The line tells us the position of the difference (chromosome and position), the genetic sequence in the reference, and the corresponding sequence in each of the three Africans. Research is ongoing to understand the full effects of these genetic differences; some cause diseases such as Tay-Sachs and Hemophilia, while others determine your blood type and eye color.
+
+Before we start processing the file, let's get a high-level view of the file that we're about to work with.
+
+What's the file size (in kilo-bytes), and how many lines are in the file?
+(<a href="javascript:toggle('q3_h1');">Click for a hint about the file size</a>)
+<div id="q3_h1" class="answer">
+There's an option to <code>ls</code> that will print the file sizes in a more human-friendly format.
+</div>
+(<a href="javascript:toggle('q3_h2');">Click for a hint about the number of lines</a>)
+<div id="q3_h2" class="answer">
+<code>$ man wc</code>
+</div>
+(<a href="javascript:toggle('q3_a');">Click for our answers</a>)
+<div id="q3_a" class="answer">
+We should get a file size around 3.6 MB with:
+<code>$ ls -lh 1000gp.vcf</code>
+Alternatively, the command <code>du</code> can be used to achieve a similar result:
+<code>$ du -h 1000gp.vcf</code>
+
+We find there are 45034 lines with:
+<code>$ wc -l 1000gp.vcf</code>
+</div>
+</li>
+
+<li>
+Because this file is so large, you're going to almost always want to pipe ("|") the result of any command to <code>less</code> (a simple text viewer, type 'q' to exit) or <code>head</code> (to print the first 10 lines) so that you don't accidentally print 45,000 lines to the screen.
+
+Let's start by printing the first 5 lines to see what it looks like.
+(<a href="javascript:toggle('q3b_a');">Click for our answer</a>)
+<div id="q3b_a" class="answer">
+<code>$ head -5 1000gp.vcf</code>
+</div>
+</li>
+
+
+<li>
+That isn't very interesting; it's just a bunch of the comments at the beginning of the file (they all start with "#")! Print the first 20 lines to see more of the file.
+(<a href="javascript:toggle('q4_a');">Click for our answer</a>)
+<div id="q4_a" class="answer">
+<code>$ head -20 1000gp.vcf</code>
+</div>
+</li>
+
+
+<li>
+Okay, so now we can see the basic structure of the file. A few comment lines that start with "#" or "##" and then a bunch of lines of data that contain all the data and are pretty hard to understand. Each line of data contains the same number of fields, and all fields are separated with TABs. These fields are:
+<ol>
+<li>the chromosome (which volume the difference is in)</li>
+<li>the position (which character in the volume the difference starts at)</li>
+<li>the ID of the difference</li>
+<li>the sequence in the reference human(s)</li>
+</ol>
+The rest of the columns tell us, in a rather complex way, a bunch of additional information about that position, including: the predicted sequence for each of the three Africans and how confident the scientists are that these sequences are correct.
+
+To start analyzing the actual data, we have to remove the header. How can we print the first 10 non-header lines (those that <em>don't</em> start with a "#")?
+(<a href="javascript:toggle('q5_h1');">Click for a hint</a>)
+<div id="q5_h1" class="answer">
+<code>$ man grep</code>
+</div>
+(<a href="javascript:toggle('q5_h2');">Click for another hint</a>)
+<div id="q5_h2" class="answer">
+You can use a pipe ("|") to connect the output of <code>grep</code> to the input of <code>head</code>.
+</div>
+(<a href="javascript:toggle('q5_a');">Click for our answer</a>)
+<div id="q5_a" class="answer">
+<code>$ grep -v "^#" 1000gp.vcf | head</code>
+
+Why are neither of these correct?
+<code>$ grep -v "#" 1000gp.vcf | head</code>
+<code>$ grep -v "^##" 1000gp.vcf | head</code>
+</div>
+</li>
+
+
+<li>
+How many lines of data are in the file (rather than counting the number of header lines and subtracting, try just counting the number of data lines)?
+(<a href="javascript:toggle('q6_h');">Click for a hint</a>)
+<div id="q6_h" class="answer">
+Instead of piping to <code>head</code>, try piping to <code>wc</code>.
+</div>
+(<a href="javascript:toggle('q6_a');">Click for our answer</a>)
+<div id="q6_a" class="answer">
+<code>$ grep -v "^#" 1000gp.vcf | wc -l</code>
+should print <code>45024</code>
+</div>
+</li>
+
+
+<li>
+Where these differences are located can be important. If all the differences between two encyclopedias were in just the first volume, that would be interesting. The first field of each data line is the name of the chromosome that the difference occurs on (which volume we're on). Print the first 10 chromosomes, one per line.
+(<a href="javascript:toggle('q7_h1');">Click for a hint</a>)
+<div id="q7_h1" class="answer">
+You can extract a column from a tab-delimited text file using the <code>cut</code> command.
+</div>
+(<a href="javascript:toggle('q7_h2');">Click for another hint</a>)
+<div id="q7_h2" class="answer">
+Use <code>grep</code> to print only non-comment lines, and <code>cut</code> to extract the chromosome column.
+</div>
+(<a href="javascript:toggle('q7_a');">Click for our answer</a>)
+<div id="q7_a" class="answer">
+<code>$ grep -v "^#" 1000gp.vcf | cut -f 1 | head</code>
+</div>
+</li>
+
+<li>
+As you should have observed, the first 10 lines are on numbered chromosomes. Every normal cell in your body has 23 pairs of chromosomes, 22 pairs of "autosomal" chromosomes (these are numbered 1-22) and a pair of sex chromosomes (two Xs if you're female, an X and a Y if you're male). If you've heard of the genetics company <a href="https://www.23andme.com">23andMe</a>, the 23 refers to these 23 pairs of chromosomes.
+
+Let's look at which chromosomes these variations are on. Print a list of the chromosomes that are in the file (each chromosome name should only be printed once, so you should only print 23 lines).
+(<a href="javascript:toggle('q8_h1');">Click for a hint</a>)
+<div id="q8_h1" class="answer">
+You need to remove all the duplicate lines from your previous answer.
+</div>
+(<a href="javascript:toggle('q8_h2');">Click for a hint</a>)
+<div id="q8_h2" class="answer">
+<code>sort</code> has an option that should make this easier.
+</div>
+(<a href="javascript:toggle('q8_a');">Click for our answer</a>)
+<div id="q8_a" class="answer">
+<code>$ grep -v "^#" 1000gp.vcf | cut -f 1 | sort -u</code>
+</div>
+</li>
+
+
+<li>
+Rather than using <code>sort</code> to print unique results, a common pipeline is to first sort and then pipe to another UNIX command, <code>uniq</code>. The <code>uniq</code> command takes <em>sorted</em> input and prints only unique lines, but it provides more flexibility than just using <code>sort</code> by itself. Keep in mind, if the input isn't sorted, <code>uniq</code> won't work properly.
+
+Using <code>sort</code> and <code>uniq</code>, print the number of times each chromosome occurs in the file.
+(<a href="javascript:toggle('q9_h1');">Click for a hint</a>)
+<div id="q9_h1" class="answer">
+<code>$ man uniq</code>
+</div>
+(<a href="javascript:toggle('q9_h2');">Click for another hint</a>)
+<div id="q9_h2" class="answer">
+Instead of using <code>sort</code> to remove duplicates, just use it to sort and pipe the result to <code>uniq</code>.
+</div>
+(<a href="javascript:toggle('q9_a');">Click for our answer</a>)
+<div id="q9_a" class="answer">
+<code>$ grep -v "^#" 1000gp.vcf | cut -f 1 | sort | uniq -c</code>
+</div>
+</li>
+
+
+<li>
+Add to your previous solution to list the chromosomes from most frequently observed to least frequently observed.
+(<a href="javascript:toggle('q10_h1');">Click for a hint</a>)
+<div id="q10_h1" class="answer">
+<code>$ man sort</code>
+</div>
+(<a href="javascript:toggle('q10_h2');">Click for another hint</a>)
+<div id="q10_h2" class="answer">
+Make sure you're sorting in descending order. By default, <code>sort</code> sorts in ascending order.
+</div>
+(<a href="javascript:toggle('q10_a');">Click for our answer</a>)
+<div id="q10_a" class="answer">
+<code>$ grep -v "^#" 1000gp.vcf | cut -f 1 | sort | uniq -c | sort -r</code>
+should output the following:
+<pre>
+3721 2
+3387 1
+3224 4
+3219 3
+2894 5
+2860 6
+2527 8
+2525 7
+2203 10
+2166 11
+2032 12
+1865 9
+1656 13
+1409 14
+1362 16
+1304 X
+1275 18
+1265 15
+1097 17
+ 993 20
+ 814 19
+ 661 21
+ 565 22
+</pre>
+
+* Because <code>uniq -c</code> right-aligns the counts, we don't have to specify the <code>-n</code> flag to <code>sort</code>, although doing so would also be correct.
+</div>
+</li>
+
+
+<li>
+The autosomal chromosomes (1-22) are named according to their size. The largest of them is chromosome 1, while the smallest is chromosome 22. Does it look like differences occur relatively randomly across the genome, or are some chromosomes more different than you'd expect at random (very roughly taking their sizes into account)?
+(<a href="javascript:toggle('q10b_h1');">Click for a hint</a>)
+<div id="q10b_h1" class="answer">
+It's worth noting that the chromosomes were numbered by the sizes of the actual molecules, not how much of them had been sequenced.
+
+Wikipedia has a nice table of chromosome sizes and how much of each has been sequenced (and you can sort it):
+<a href="http://en.wikipedia.org/wiki/Human_chromosome#Human_chromosomes">http://en.wikipedia.org/wiki/Human_chromosome#Human_chromosomes</a>
+
+Notice anything?
+</div>
+(<a href="javascript:toggle('q10b_a');">Click for our hypothesis</a>)
+<div id="q10b_a" class="answer">
+Since variation can only be found in the known sequence, the order you printed corresponds closely to ordering by the number of bases sequenced (rather than the total number of bases).
+
+Given this, it seems like differences occur relatively randomly across the genome. We see more differences on longer chromosomes, fewer on shorter, without any striking outliers.
+</div>
+</li>
+
+<li>
+This is great, but biologists might also like to see the chromosomes ordered by their number (not dictionary order), since different chromosomes have different attributes and this ordering allows them to find a specific chromosome more easily.
+(<a href="javascript:toggle('q11_h');">Click for another hint</a>)
+<div id="q11_h" class="answer">
+A lot of the power of <code>sort</code> comes from the fact that you can specify which fields to sort on, and the order in which to sort them. In this case you only need to sort on one field.
+</div>
+(<a href="javascript:toggle('q11_a');">Click for our answer</a>)
+<div id="q11_a" class="answer">
+<code>$ grep -v "^#" 1000gp.vcf | cut -f 1 | sort | uniq -c | sort -k 2n</code>
+</div>
+</li>
+</ul>
+
+
+<h2>Exercise Part 3 (scripts and svn)</h2>
+<ul>
+<li>
+Wonderful! Now we have a (long) command for printing chromosome statistics from our <code>1000gp.vcf</code> file. Using <code>nano</code>, create a new file, "chrom_stats.sh", with just your answer to the previous question in it.
+(<a href="javascript:toggle('q12a_a');">Click for our answer</a>)
+<div id="q12a_a" class="answer">
+Type the following to open a new file:
+<code>$ nano chrom_stats.sh</code>
+Type in the command. Type ^o to save and ^x (where ^ means the control key).
+</div>
+</li>
+
+<li>
+Just to be illustrate the flexibility of the shell, try creating the same file directly from the shell (without a text editor). Once you do, you can use <code>cat</code> to make sure the contents of the file are exactly what you expect.
+(<a href="javascript:toggle('q12b_h1');">Click for a hint</a>)
+<div id="q12b_h1" class="answer">
+You can use <code>echo</code> to print something and <code>></code> to redirect to a file.
+</div>
+(<a href="javascript:toggle('q12b_h2');">Click for another hint</a>)
+<div id="q12b_h2" class="answer">
+Since our long command has double-quotes in it, you either need to use single-quotes or escape these with back-slashes.
+</div>
+(<a href="javascript:toggle('q12b_a');">Click for our answer</a>)
+<div id="q12b_a" class="answer">
+<code>$ echo 'grep -v "^#" 1000gp.vcf | cut -f 1 | sort | uniq -c | sort -k 2n' > chrom_stats.sh</code>
+<code>$ cat chrom_stats.sh</code>
+</div>
+</li>
+
+
+<li>
+Now, execute your new script to print the chromosome statistics.
+(<a href="javascript:toggle('q13_h1');">Click for a hint</a>)
+<div id="q13_h1" class="answer">
+You may have to change the permissions to allow you to execute it.
+</div>
+(<a href="javascript:toggle('q13_h2');">Click for a hint</a>)
+<div id="q13_h2" class="answer">
+It's good form to only make permissions as permissive as necessary. So, rather than allow everyone to execute the file, it is better to just allow you to execute it.
+</div>
+(<a href="javascript:toggle('q13_a');">Click for our answer</a>)
+<div id="q13_a" class="answer">
+<code>$ chmod u+x chrom_stats.sh
+$ ./chrom_stats.sh
+</code>
+
+Note that it is <code>u+x</code> instead of just <code>+x</code> or <code>a+x</code>. This only adds the ability for the owner to execute it, whereas the other two options would allow anyone to execute it.
+</div>
+</li>
+
+
+<li>
+We'd like to be able to use this script in the future with arbitrary VCF files, instead of just our <code>1000gp.vcf</code> file. Edit the script so that it takes VCF-formatted text input on stdin and prints out chromosome statistics on stdout. This is simpler than you might think.
+(<a href="javascript:toggle('q14_h');">Click for a hint</a>)
+<div id="q14_h" class="answer">
+If <code>grep</code> isn't given an input file, it will read from stdin.
+</div>
+(<a href="javascript:toggle('q14_a');">Click for our answer</a>)
+<div id="q14_a" class="answer">
+Change
+<code>grep -v "^#" 1000gp.vcf | ...</code>
+to
+<code>grep -v "^#" | ...</code>
+
+* Since this is in a file instead of the shell prompt, we aren't showing the "$" at the beginning of the line.
+</div>
+</li>
+
+<li>
+Now that we have a script that reads from stdin and prints to stdout, how do we run it on the <code>1000gp.vcf</code> file to get the same output as before?
+(<a href="javascript:toggle('q15_h1');">Click for a hint</a>)
+<div id="q15_h1" class="answer">
+The <code>cat</code> command is used to print files to stdout.
+</div>
+(<a href="javascript:toggle('q15_h2');">Click for another hint</a>)
+<div id="q15_h2" class="answer">
+You can pipe the output of <code>cat</code> directly into our script.
+</div>
+(<a href="javascript:toggle('q15_h3');">Click for yet another hint</a>)
+<div id="q15_h3" class="answer">
+Just like before, in order to tell the shell that the <code>chrom_stats.sh</code> file we want to execute is the one in our current directory, we need to use <code>./chrom_stats.sh</code>.
+</div>
+(<a href="javascript:toggle('q15_a');">Click for our answer</a>)
+<div id="q15_a" class="answer">
+<code>$ cat 1000gp.vcf | ./chrom_stats.sh</code>
+</div>
+</li>
+
+<li>
+Finally, add a copy of this file to your folder in the class SVN repository.
+(<a href="javascript:toggle('q16_a');">Click for our answer</a>)
+<div id="q16_a" class="answer">
+<ol>
+<li><code>cp chrom_stats.sh /path/to/repo/participants/user/</code></li>
+<li>Add the file to subversion version control</li>
+<li>Commit your changes</li>
+</ol>
+</div>
+</li>
+
+<li>
+<b>Fin.</b> Comments, questions, and suggestions are encouraged and appreciated.
+Thanks to Tommy Guy, Jon Pipitone, Greg Wilson, and Elango Cheran for their help with these exercises.
+</li>
\ No newline at end of file
projects / swc-workshop.git / commitdiff