[Jon_Keats]

The Shell Script from Hell or My Masterpiece you Decide

I'll preface this with acknowledging all the help from various people who have replied to this thread and I'll promise to try and figure out perl over shell scripts, but two weeks ago the "Hello World!" script in any programing language was "Greek" to me so bare with my ignorance for a week or two more. But in the interest of sharing here is my current bwa samse pipeline. I'm stuck at not being able to merge the sorted bam files but at least this lets me get back to the lab and opens up some of my time to write the back log of papers I need to polish off so I can get my own lab and make someone else deal with these issues
**Currently, this uses all default settings so modify commands as you see fit***

#!/bin/sh

# Batch_BWAsamse_V1.sh
# Created by Jonathan Keats on 04/05/10.
# This script is designed to take a batch of raw Illumina 1.3+ reads to sorted BAM files.
# It is designed to be initiated from a folder called "NGS" with a specific subdirectory structure
# Use the script called "Create_NGS_DirectoryStructureV1.sh" or the other variations (see Jon_Keats SeqAnswers Introduction Thread) to create the analysis directory structure
# To run this script you MUST first place your reference file in "NGS/RefGenomes/BWA_Indexed" and have run the "bwa index" command to create the BWT index files
# If you are NOT using a reference genome called "hg18.fa" you will NEED to edit lines 124 and 150 accordingly!
# The script is based on having ***RENAMED*** Illumina files in NGS/InputSequence/Illumina/Read1 and if available NGS/InputSequence/Illumina/Read2 if paired end runs
# The renamed format MUST be "YourSampleName_RunX_LaneX_R1.txt" and "YourSampleName_RunX_LaneX_R2.txt" otherwise files will be overwritten, paired-end analysis will not be possible, and lane merging will not be possible 
# At each step it queries specific folders for available files and passes them to the next analysis module
# After each step the filename extension of the output files are corrected. (ie. "MySequenceFile.txt.fastq" to "MySequenceFile.fastq")
# Order of Embedded Steps	- Converts Illumina 1.3+ fastq files "s_1_sequence.txt" to Sanger fastq using "maq ill2sanger" command
#							- Aligns created fastq files to reference genome using "bwa aln" command
#							- Generates SAM files from alignment files using "bwa samse" command
#							- Converts SAM files to BAM files using "samtools view" command
#							- Sorts BAM files using "samtools sort" command
#							- Final output files are archived then the input and analysis directories are cleaned-up and readied for the next analysis batch
#    ***DURING CLEAN-UP THE ENTIRE CONTENTS OF SEVERAL FOLDERS ARE DELETED, SO ONLY PLACE ADDITIONAL FILES IN THE "NGS/FINALOUTPUTS" DIRECTORY OR OUTSIDE THE NGS FOLDER SUBDIRECTORIES***
# The script creates a log file in "NGS/AnalysisNotes" to track the steps completed and the time each step started and finished
# Some of the log events will print to both the terminal screen and the log file so you can see what is going on
# On our workstation this will process two lanes of PE data overnight or a full run over the weekend

# Sould you merge files (lane 1 and 2) before analysis as illumina fastq files or after sorting BAM files?
# If you merge before analysis does this mess up the a paired-end analysis?

# Carabelle - Can you have a fail safe start query to kill the script if it is not initiated from the "NGS" folder? Probably doesn't run anyways
# Carabelle - Also can you check if destination folders are empty and if not kill the script?
# Carabelle - Such as if directory NGS/SAMfiles is NOT empty kill script and print "The NGS/SAMfiles directory is not empty please remove files otherwise they will be deleted"
# Ligne=Line en francais (I had a significant amount of help from our French post-docs husband who is a unix programer in Paris)

#Starting directory = NGS/

echo ***Starting Analysis Batch***
date '+%m/%d/%y %H:%M:%S'

#The following step creates the log file in the AnalysisNotes subdirectory the first time the script is run
#On subsequent runs the results are printed at the bottom of the pre-existing log file

echo ***Starting Analysis Batch*** >> AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> AnalysisNotes/Analysis.log

#In this step we will converted the original Illumina 1.3+ fastq format files to Sanger fastq format
#Input files from NGS/InputSequences/Illumina/Read1 and NGS/InputSequences/Illumina/Read2
#Output files to NGS/InputSequences/SangerFastq

echo Starting Step1a - Read1 Fastq Conversion with maq ill2sanger
date '+%m/%d/%y %H:%M:%S'
echo Starting Step1a - Read1 Fastq Conversion with maq ill2sanger >> AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> AnalysisNotes/Analysis.log
cd InputSequence/Illumina/Read1
#Current directory = NGS/InputSequence/Illumina/Read1
echo Converting the following Illumina files from the Read1 folder:
for ligne in `ls *.txt`
do                                                                     
echo $ligne
done
echo Converting the following Illumina files from the Read1 folder: >> ../../../AnalysisNotes/Analysis.log
for ligne in `ls *.txt`
do
echo $ligne >> ../../../AnalysisNotes/Analysis.log
done
for ligne in `ls *.txt`
do
maq ill2sanger $ligne ../../SangerFastq/$ligne.fastq
done
echo Finished Step1a - Read1 Fastq Conversion
date '+%m/%d/%y %H:%M:%S'
echo Finished Step1a - Read1 Fastq Conversion >> ../../../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../../../AnalysisNotes/Analysis.log
cd ../Read2
#Current directory = NGS/InputSequence/Illumina/Read2
echo Starting Step1b - Read2 Fastq Conversion with maq ill2sanger
date '+%m/%d/%y %H:%M:%S'
echo Starting Step1b - Read2 Fastq Conversion with maq ill2sanger >> ../../../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../../../AnalysisNotes/Analysis.log
echo Converting the following Illumina files from the Read2 folder:
for ligne in `ls *.txt`
do                                                                     
echo $ligne
done
echo Converting the following Illumina files from the Read2 folder: >> ../../../AnalysisNotes/Analysis.log
for ligne in `ls *.txt`
do
echo $ligne >> ../../../AnalysisNotes/Analysis.log
done
for ligne in `ls *.txt`
do
maq ill2sanger $ligne ../../SangerFastq/$ligne.fastq
done
cd ../../SangerFastq/
#Current directory = NGS/InputSequence/SangerFastq
#The next step will change the file extensions of ths Sanger Fastq outputs from "File_Run1_Lane1_R1.txt.fastq" to "File_Run1_Lane1_R1.fastq"
old_ext=txt.fastq
new_ext=fastq
find . -type f -name "*.$old_ext" -print | while read file
do
    mv $file ${file%${old_ext}}${new_ext}
done
echo Finished Step1b - Read2 Fastq Conversion
date '+%m/%d/%y %H:%M:%S'
echo Finished Step1b - Read2 Fastq Conversion >> ../../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../../AnalysisNotes/Analysis.log

#In the next step we will align the converted sanger fastq format files to the reference genome (hg18.fa)

echo Starting Step2 - bwa aln process
date '+%m/%d/%y %H:%M:%S'
echo Starting Step2 - bwa aln process >> ../../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../../AnalysisNotes/Analysis.log
echo The following fastq files will be aligned:
for ligne in `ls *.fastq`
do                                                                     
echo $ligne
done
echo The following fastq files will be aligned: >> ../../AnalysisNotes/Analysis.log
for ligne in `ls *.fastq`
do
echo $ligne >> ../../AnalysisNotes/Analysis.log
done
for ligne in `ls *.fastq`
do
bwa aln ../../RefGenomes/BWA_Indexed/hg18.fa $ligne > $ligne.sai 	 
done
echo Finished Step2 - bwa aln process
date '+%m/%d/%y %H:%M:%S'
echo Finished Step2 - bwa aln process >> ../../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../../AnalysisNotes/Analysis.log

#In the next step we will generate SAM files for the alignments using bwa samse

echo Starting Step3 - bwa samse process
date '+%m/%d/%y %H:%M:%S'
echo Starting Step3 - bwa samse process >> ../../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../../AnalysisNotes/Analysis.log
echo The following alignment files will be converted to SAM files:
for ligne in `ls *.sai`
do                                                                     
echo $ligne
done
echo The following alignment files will be converted to SAM files: >> ../../AnalysisNotes/Analysis.log
for ligne in `ls *.sai`
do
echo $ligne >> ../../AnalysisNotes/Analysis.log
done
#To get this step to work I had to trick the system to get the conventional (bwa samse <database.fasta> <aln.sai> <input.fastq> > aln.sam) command as I can only pull up one list at a time to automatically feed files into the command
for ligne in `ls *.fastq`
do
bwa samse ../../RefGenomes/BWA_Indexed/hg18.fa $ligne.sai $ligne > ../../SAMfiles/$ligne.sam
done
#This step renames the bwa aln output files from "File_Run1_Lane1_R1.fastq.sai" to "File_Run1_Lane1_R1.sai"
#This is done on purpuse at this point to allow the preceeding step to work
old_ext=fastq.sai
new_ext=sai
find . -type f -name "*.$old_ext" -print | while read file
do
    mv $file ${file%${old_ext}}${new_ext}
done 
cd ../../SAMfiles/
#Current directory = NGS/SAMfiles
#This step renames the bwa samse output files from "File_Run1_Lane1_R1.fastq.sam" to "File_Run1_Lane1_R1_bwase.sam"
#We up in the "_bwase" to indicate the files are bwa alignments converted to SAM files with the bwa samse command
old_ext=.fastq.sam
new_ext=_bwase.sam
find . -type f -name "*$old_ext" -print | while read file
do
    mv $file ${file%${old_ext}}${new_ext}
done 
echo Finished Step3 - bwa samse process
date '+%m/%d/%y %H:%M:%S'
echo Finished Step3 - bwa samse process >> ../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../AnalysisNotes/Analysis.log

#In the next step we will convert each SAM file to a BAM file

echo Starting Step4 - samtools SAM to BAM conversion
date '+%m/%d/%y %H:%M:%S'
echo Starting Step4 - samtools SAM to BAM conversion >> ../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../AnalysisNotes/Analysis.log
echo The following SAM files will be converted to BAM files:
for ligne in `ls *.sam`
do                                                                     
echo $ligne
done
echo The following SAM files will be converted to BAM files: >> ../AnalysisNotes/Analysis.log
for ligne in `ls *.sam`
do                                                                     
echo $ligne >> ../AnalysisNotes/Analysis.log
done
for ligne in `ls *.sam`
do
samtools view -bS -o ../BAMfiles/Original/$ligne.bam $ligne
done
cd ../BAMfiles/Original
#Current directory = NGS/BAMfiles/Original
#This step renames the samtools view output files from "File_Run1_Lane1_R1_bwase.sam.bam" to "File_Run1_Lane1_R1_bwase.bam"
old_ext=sam.bam
new_ext=bam
find . -type f -name "*.$old_ext" -print | while read file
do
    mv $file ${file%${old_ext}}${new_ext}
done 
echo Finished Step4 - samtools SAM to BAM conversion
date '+%m/%d/%y %H:%M:%S'
echo Finished Step4 - samtools SAM to BAM conversion >> ../../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../../AnalysisNotes/Analysis.log

#In the next step we will sort the BAM file by chromosome coordinate

echo Starting Step4 - samtools BAM sorting process
date '+%m/%d/%y %H:%M:%S'
echo Starting Step4 - samtools BAM sorting process >> ../../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../../AnalysisNotes/Analysis.log
echo The following BAM files will be sorted:
for ligne in `ls *.bam`
do                                                                     
echo $ligne
done
echo The following BAM files will be sorted: >> ../../AnalysisNotes/Analysis.log
for ligne in `ls *.bam`
do                                                                     
echo $ligne >> ../../AnalysisNotes/Analysis.log
done
for ligne in `ls *.bam`
do                                                                     
samtools sort $ligne ../Sorted/$ligne
done
cd ../Sorted
#Current directory = NGS/BAMfiles/Sorted
#This step renames the samtools sort output files from "File_Run1_Lane1_R1_bwase.bam.bam" to "File_Run1_Lane1_R1_bwase_sorted.bam
old_ext=.bam.bam
new_ext=_sorted.bam
find . -type f -name "*$old_ext" -print | while read file
do
    mv $file ${file%${old_ext}}${new_ext}
done 
echo Finished Step4 - samtools BAM sorting process
date '+%m/%d/%y %H:%M:%S'
echo Finished Step4 - samtools BAM sorting process >> ../../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../../AnalysisNotes/Analysis.log
cd ../..
#Current directory = NGS/

#In the next step we will move files from the input and analysis directories to their respective "NGS/FinalOutputs" subdirectory so the pipeline can be used another time.

echo Cleaning up Input and Analysis Directories
date '+%m/%d/%y %H:%M:%S'
echo Cleaning up Input and Analysis Directories >> AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> AnalysisNotes/Analysis.log
cd InputSequence/Illumina/Read1
#Current directory = NGS/InputSequence/Illumina/Read1
echo Moving the following Illumina Fastq Files from NGS/InputSequence/Illumina/Read1 to NGS/FinalOutputs/Illumina/Read1:
for ligne in `ls *.txt`
do
echo $ligne
done
echo Moving Illumina Fastq Files from NGS/InputSequence/Illumina/Read1 to NGS/FinalOutputs/Illumina/Read1 >> ../../../AnalysisNotes/Analysis.log
for ligne in `ls *.txt`
do
echo $ligne >> ../../../AnalysisNotes/Analysis.log
done
for ligne in `ls *.txt`
do
mv $ligne ../../../FinalOutputs/Illumina/Read1/
done
cd ../Read2
#Current directory = NGS/InputSequence/Illumina/Read2
echo Moving the following Illumina Fastq Files from NGS/InputSequence/Illumina/Read2 to NGS/FinalOutputs/Illumina/Read2:
for ligne in `ls *.txt`
do
echo $ligne
done
echo Moving Illumina Fastq Files from NGS/InputSequence/Illumina/Read2 to NGS/FinalOutputs/Illumina/Read2 >> ../../../AnalysisNotes/Analysis.log
for ligne in `ls *.txt`
do
echo $ligne >> ../../../AnalysisNotes/Analysis.log
done
for ligne in `ls *.txt`
do
mv $ligne ../../../FinalOutputs/Illumina/Read2/
done
echo Moving Illumina Input Files Complete
date '+%m/%d/%y %H:%M:%S'
echo Moving Illumina Input Files Complete >> ../../../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../../../AnalysisNotes/Analysis.log
cd ../../SangerFastq
#Current directory = NGS/InputSequence/SangerFastq
echo Moving the following Sanger Format Fastq Files from NGS/InputSequence/SangerFastq to NGS/FinalOutputs/SangerFastq:
for ligne in `ls *.fastq`
do
echo $ligne
done
echo Moving the following Sanger Format Fastq Files from NGS/InputSequence/SangerFastq to NGS/FinalOutputs/SangerFastq: >> ../../AnalysisNotes/Analysis.log
for ligne in `ls *.fastq`
do
echo $ligne >> ../../AnalysisNotes/Analysis.log
done
for ligne in `ls *.fastq`
do
mv $ligne ../../FinalOutputs/SangerFastq/
done
echo Moving Sanger Format Fastq Files Complete
date '+%m/%d/%y %H:%M:%S'
echo Moving Sanger Format Fastq Files Complete >> ../../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../../AnalysisNotes/Analysis.log
echo Moving the following Alignment .sai Files from NGS/InputSequence/SangerFastq to NGS/FinalOutputs/AlignmentResults:
for ligne in `ls *.sai`
do
echo $ligne
done
echo Moving the following Alignment .sai Files from NGS/InputSequence/SangerFastq to NGS/FinalOutputs/AlignmentResults: >> ../../AnalysisNotes/Analysis.log
for ligne in `ls *.sai`
do
echo $ligne >> ../../AnalysisNotes/Analysis.log
done
for ligne in `ls *.sai`
do
mv $ligne ../../FinalOutputs/AlignmentResults/
done
echo Moving Alignment Results Files Complete
date '+%m/%d/%y %H:%M:%S'
echo Moving Alignment Results Files Complete >> ../../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../../AnalysisNotes/Analysis.log
cd ../../BAMfiles/Sorted/
#Current directory = NGS/BAMfiles/Sorted/
echo Moving the following Sorted BAM files from NGS/BAMfiles/Sorted/ to NGS/FinalOutputs/SortedBAMfiles:
for ligne in `ls *.bam`
do
echo $ligne
done
echo Moving the following Sorted BAM files from NGS/BAMfiles/Sorted/ to NGS/FinalOutputs/SortedBAMfiles: >> ../../AnalysisNotes/Analysis.log
for ligne in `ls *.bam`
do
echo $ligne >> ../../AnalysisNotes/Analysis.log
done
for ligne in `ls *.bam`
do
mv $ligne ../../FinalOutputs/SortedBAMfiles/
done
echo Moving Sorted BAM Files Complete
date '+%m/%d/%y %H:%M:%S'
echo Moving Alignment Results Files Complete >> ../../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../../AnalysisNotes/Analysis.log
cd ../../Samfiles
#Current directory = NGS/SAMfiles
echo Deleting the following SAM Files from NGS/SAMfiles:
for ligne in `ls *.sam`
do
echo $ligne
done
echo Deleting the following SAM Files from NGS/SAMfiles: >> ../AnalysisNotes/Analysis.log
for ligne in `ls *.sam`
do
echo $ligne >> ../AnalysisNotes/Analysis.log
done
for ligne in `ls *.sam`
do
rm $ligne
done
echo Deleting SAM Files Complete
date '+%m/%d/%y %H:%M:%S'
echo Deleting SAM Files Complete >> ../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../AnalysisNotes/Analysis.log
cd ../BAMfiles/Original/
#Current directory = NGS/BAMfiles/Original
echo Deleting the following BAM Files from NGS/BAMfiles/Original:
for ligne in `ls *.bam`
do
echo $ligne
done
echo Deleting the following BAM Files from NGS/BAMfiles/Original: >> ../../AnalysisNotes/Analysis.log
for ligne in `ls *.bam`
do
echo $ligne >> ../../AnalysisNotes/Analysis.log
done
for ligne in `ls *.bam`
do
rm $ligne
done
echo Deleting BAM Files from NGS/BAMfiles/Original Complete
date '+%m/%d/%y %H:%M:%S'
echo Deleting BAM Files from NGS/BAMfiles/Original Complete >> ../../AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> ../../AnalysisNotes/Analysis.log
cd ../..
#Current directory = NGS/
echo Cleaning up Input and Analysis Directories Complete
date '+%m/%d/%y %H:%M:%S'
echo Cleaning up Input and Analysis Directories Complete >> AnalysisNotes/Analysis.log
date '+%m/%d/%y %H:%M:%S' >> AnalysisNotes/Analysis.log
echo ***Analysis Batch Complete***
echo ***Analysis Batch Complete*** >> AnalysisNotes/Analysis.log

# Future development: Merge BAM files based on filename structure "RunX_LaneX_RX", Index files for IGV viewing
# Figure out Perl or Python so I don't look like such an amateur :)

[Jon_Keats]

To Genome or Not to Genome

I'm sure I'm not the only person facing a conundrum on which genome version to use for alignment. In general I'm a believer in the "out with the old, in with the new" mentality so my gut is pushing me to use the newest genome build GRCh37/hg19/Ensembl55+ over the old version NCBI36/hg18/Ensembl38-54. Plus CCDS has been on GRCh37 for nearly a year now and they just released dbSNP131 on GRCh37 so it's getting harder and harder to ignore the newest build. BUT, for my dataset I have paired copy number data on Agilent arrays, which still use hg18 design files. So after wasting a couple of days aligning things to GRCh37 (1000 genomes version) I'm going back to redo it with hg18 so my coordinates match between the two datasets, though wasting might be the wrong word as I'm betting I'll want to be on the new version in a month. To get this genome reference file I decided to use the UCSC chromosome files (http://hgdownload.cse.ucsc.edu/golde...8/chromosomes/) since I'd seen someone mentioning the advantage of their simple naming structure. I did not use the various haplotype builds or the random contigs.
- Downloaded each individual chromosome file Chr1-Chr22, ChrX, ChrY, and ChrM.
- Placed them all in my "ngs/refgenomes/genome_downloads/hg18" folder
- Navigate to the folder in Terminal "cd ngs/refgenomes/genome_downloads/hg18"
- Decompress the files using gunzip, type "gunzip *.fa.gz"
- Then merge the files using unix cat command [cat (infile1) (infile2) (infileN) > (outfile)]
- Type "cat chr1.fa chr2.fa chr3.fa chr4.fa chr5.fa chr6.fa chr7.fa chr8.fa chr9.fa chr10.fa chr11.fa chr12.fa chr13.fa chr14.fa chr15.fa chr16.fa chr17.fa chr18.fa chr19.fa chr20.fa chr21.fa chr22.fa chrX.fa chrY.fa chrM.fa > hg18.fasta"
- check the ouput file contains all chromosomes in the correct order using unix grep command to pull out the fasta header lines
- Type "grep ">" hg18.fa"
- Which should output the following if correct:
>chr1
>chr2
>chr3
>chr4
>chr5
>chr6
>chr7
>chr8
>chr9
>chr10
>chr11
>chr12
>chr13
>chr14
>chr15
>chr16
>chr17
>chr18
>chr19
>chr20
>chr21
>chr22
>chrX
>chrY
>chrM

A similar process should be possible with any genome version or source like ensembl. Just make sure the end platform supports the version you are using. And be careful with the annotation files if you get them from UCSC as the positions are 0 base indexed.

[Bruins]

Hi all,

Let me do my best to help out with the Perl.

A lot of basic operations have already been written in Perl. You can easily use these bits of code. If you need something done, check http://search.cpan.org/. This goes for bioinformatics things as well: check http://www.bioperl.org for stuff like dealing with sequence files.

For your Perl script I used File::Path (http://search.cpan.org/~jhi/perl-5.8.0/lib/File/Path.pm). It took me about 5 minutes to write.

#!/usr/bin/perl
use strict;
use File::Path; # CPAN modules are this easy to use

# the dirs
my $dirs = [ 'NGS/AnalysisNotes',
	  'NGS/ApplicationDownloads',
	  'NGS/InputSequence',
	  'NGS/SAMfiles', 
	  'NGS/BAMfiles/Merged', 'NGS/BAMfiles/Original', 'NGS/BAMfiles/Sorted', 
	  'NGS/FinalOutputs/AlignmentResults', 'NGS/FinalOutputs/Illumina', 'NGS/FinalOutputs/SangerFasta', 'NGS/FinalOutputs/SortedBAMfiles', 'NGS/FinalOutputs/MergedBAMfiles', 
	  'NGS/RefGenomes/BFAST_Indexed', 'NGS/RefGenomes/BOWTIE_Indexed', 'NGS/RefGenomes/GenomeDownloads',
	  'NGS/Scripts/ScriptBackups' ];

print "***Creating Pipeline Directory Structure***\n";

# eval will evaluate whatever is inside the {}. Any errors are stored in $@
# mkpath is the File::Path method that will create the given directories
# it takes three arguments in total: the dirs, whether you wish to print the created directories (in that case, 1) or not (0) and the last argument has to do with file permissions, I omitted it
eval { mkpath ( $dirs, 1 ); };

# now check if there were errors. If so, print.
if ( $@ ) {
  print "Couldn't create NGS directories: $@\n";
}
print "Pipeline Directory Structure Created\n";

Hope this helps?

Cheers,
Wil

[krobison]

Better metrics have to do with readability & maintainability rather than utter brevity -- it is quite possible to write extremely short code in Perl which looks like hieroglyphics.

Your program is straightforward -- but as it grows it will get less so. One thing in particular I addressed is if you change the name of certain of your directories -- in your shell script some of the directories are listed twice, which is an opportunity for error should you decide to rename them.

P.S. Only a terminal -- no internal -- Asparagine on my name!

[Jon_Keats]

Okay still working on figuring out the Perl version but it seem pretty obvious less some of the perl coding language. However, I'm nothing if not competitive so here is the three line version. Not sure if "krobinson" perl version is fewer characters or not but this must win for number of lines.
I guess reading the manuals for functions helps mkdir -p sure helps! I'll work on figuring out git version control and a bit more perl but I can't wait to see krobinson's response to the 400+ pipeline shell script. At least I took the suggestion to ditch the upper and lower case filenames, so everything is now lower case.

NOTE - THIS SCRIPT HAS BEEN UPATED TO VERSION 3 IN A SUBSEQUENT POST

#!/bin/sh

#Create_NGS_DirectoryStructureV2.sh
#Written by Jonathan Keats 04/12/2010

echo ***Creating Pipeline Directory Structure***
mkdir -p ngs/analysisnotes ngs/applicationdownloads ngs/bamfiles/merged ngs/bamfiles/original ngs/bamfiles/sorted ngs/finaloutputs/alignmentresults ngs/finaloutputs/illumina/read1 ngs/finaloutputs/illumina/read2 ngs/finaloutputs/mergedbamfiles ngs/finaloutputs/sangerfastq ngs/finaloutputs/sortedbamfiles ngs/inputsequences/illumina/read1 ngs/inputsequences/illumina/read2 ngs/inputsequences/sangerfastq ngs/refgenomes/bfast_indexed ngs/refgenomes/bowtie_indexed ngs/refgenomes/bwa_indexed ngs/refgenomes/genomedownloads ngs/samfiles ngs/scripts/scriptbackups
cp Create_NGS_DirectoryStructureV2.sh ngs/scripts/scriptbackups/
echo ***Pipeline Directory Structure Created***

[krobison]

Aaah, I'm too big a softie. Also, IMHO the perl is more readable than shell to understand what is going on (yikes! something less readable than perl!)

Personally, I prefer directory names in all lowercase. Brevity is valuable too -- "bam" & "sam" would be my choices for those directories & I probably wouldn't bother having separate ones anyway. But that is all personal preference & you should pick what suits you.

Also: DO NOT try to use a subdirectory to backup your scripts -- it is easy to install & learn git (or a similar system) to control versions. It is much too easy to accidentally copy your file the wrong way and wipe something out; version control programs track all your changes & allow you to safely move backwards & forwards through the history of your code. They can also be a handy way to replicate code across machines (and keep them in sync) as well as a path towards collaborative coding.

#!/usr/bin/perl
use strict;
print "***Creating Pipeline Directory Structure***\n"; # \n is newline
my $root="NGS";
system("mkdir $root");

# this line only makes directories with no subdirectories
&makeSubdirs($root,
                          "AnalysisNotes", "ApplicationDownloads",
                          "InputSequence","SAMfiles");
&makeSubdirs("$root/BAMfiles","Merged","Original","Sorted");
&makeSubdirs("$root/FinalOutputs",
                     "AlignmentResults","Illumina","SangerFasta",
                       "SortedBAMfiles","MergedBAMfiles");
&makeSubdirs("$root/RefGenomes",
                   "BFAST_Indexed","BOWTIE_Indexed","GenomeDownloads");
&makeSubdirs("$root/Scripts","ScriptBackups");  # DANGEROUS!!!
print "Pipeline Directory Structure Created\n";

sub makeSubdirs
{
  my @subdirs=@_;
  my $parentDir=shift(@subdirs);
  mkdir($parentDir) if (! -d $parentDir);
  foreach my $subDir(@subdirs)
  {
     mkdir("$parentDir/$subDir");
  }
}

#cd ../Scripts
#mkdir ScriptBackups
#cd ../..
#cp Create_NGS_DirectoryStructureV1.sh NGS/Scripts/ScriptBackups/
#cd NGS/
#pwd

[Jon_Keats]

Sounds like I should finish working through the "Unix and Perl for biologist" course to see if I can figure Perl out. I managed to figure out a solution for my problem about half way through the Unix section so I started with that but maybe I shouldn't have jumped the gun so much. I'll see what I can come up with this weekend.

[krobison]

Do yourself a favor & swear off shell scripts for Perl/Python/Ruby -- once you start trying to have any sort of conditional programming, shell becomes a nightmare.

Converting your script above to Perl is pretty easy (but left as an exercise for the student :-)

[Kurt]

I'm LMAO since it appears that I'm coming to these revelations at just about the exact same time you were/are. It's definitely comforting to see somebody is figuring out the exact same stuff that I am at basically the same time I am. Everytime I talk about this stuff here, people look at me like I'm ALF or something.

[Michael.James.Clark]

I went from being a pipette jockey who did qPCR for a living to writing an algorithm for SV detection in SOLiD data and publishing a whole genome sequence in a major journal.

YOU CAN DO IT TOO!

And here's how:

[Jon_Keats]

The following script will create all the directories noted in the previous post if you want to replicate the pipeline I'm putting together...

NOTE - THIS SCRIPT HAS BEEN UPDATED TO VERSION 3 IN A LATER POST

#!/bin/sh

# Create_NGS_DirectoryStructureV1.sh
# Created by Jonathan Keats on 4/5/10.
# This file will create the directory structure needed for the subsequent pipeline
# To get this script working do one of the following:
# Option 1 - Open Terminal
#		Navigate to directory of interest, "Documents" in my case (cd Documents/)
#		Type "nano Create_NGS_DirectoryStructureV1.sh"
#                        (This open the unix nano text editor)
#		Paste from "#!/bin/sh" to "echo Pipeline Directory Structure Created"
#		Control-O to save, Control-X to exit
#		Make executable, Type "chmod 755 Create_NGS_DirectoryStructureV1.sh"
# Option 2 -  Open Xcode
#		Click "File" and select "New File"
#		In "Choose a template for your new file" select "Shell Script"
#		In new file dialogue enter File Name:"Create_NGS_DirectoryStructureV1.sh"
#		Change location to directory of interest, "Documents" in my case
#		Click "Finish"
#		Paste from "#!/bin/sh" to "echo Pipeline Directory Structure Created"
#		Click "File" and select "Save"
#		Close the file, which should already be executable
# Regardless of Option used - type "./Create_NGS_DirectoryStructureV1.sh" to launch

echo ***Creating Pipeline Directory Structure***
pwd
ls
mkdir NGS
cd NGS/
mkdir AnalysisNotes
mkdir ApplicationDownloads
mkdir BAMfiles
mkdir FinalOutputs
mkdir InputSequence
mkdir RefGenomes
mkdir SAMfiles
mkdir Scripts
cd BAMfiles/
mkdir Merged
mkdir Original
mkdir Sorted
cd ../FinalOutputs/
mkdir AlignmentResults
mkdir Illumina
mkdir SangerFastq
mkdir SortedBAMfiles
mkdir MergedBAMfiles
cd Illumina/
mkdir Read1
mkdir Read2
cd ../../InputSequence/
mkdir Illumina
mkdir SangerFastq
cd Illumina/
mkdir Read1
mkdir Read2
cd ../../RefGenomes
mkdir BFAST_Indexed
mkdir BOWTIE_Indexed
mkdir BWA_Indexed
mkdir GenomeDownloads
cd ../Scripts
mkdir ScriptBackups
cd ../..
cp Create_NGS_DirectoryStructureV1.sh NGS/Scripts/ScriptBackups/
cd NGS/
pwd
ls
echo Pipeline Directory Structure Created

[Jon_Keats]

Getting setup and compiling the applications

As promised here is my next installment on getting a working environment going or at least my poor excuse of one. The first step I setup was a series of folders to manage the data off my terabyte drives and move it around as each step is completed. To make my examples more clear I've setup the following folders and subfolders:

Main working directory called "ngs" in my $HOME directory (Users/MeOrYou/) from which all steps and scripts will be called.
With primary subfolders: /ngs/analysisnotes
/ngs/applications
/ngs/bwase
/ngs/bwape
/ngs/finaloutputs
/ngs/scripts
With a number of secondary subfolders in each primary directory (See create_ngs_directorystructure_v3.sh script for full details)

I'm slowly building pipeline scripts to feed data from the input folders to final outputs that I'll try and post when complete.

The basic idea is that I have some raw reads from our Illumina GAIIs (exon capture and RNAseq PE data with each sample on two flowcell lanes) and I want to process them with BWA and view the alignments in the IGV browser. So I need to do a couple of things to the best of my understanding. Step one is to convert the Illumina raw data files, should look like "s_1_sequences.txt", to sanger fastq format files. To understand the differences please see the following references (http://en.wikipedia.org/wiki/FASTQ_format) or (http://maq.sourceforge.net/fastq.shtml) or (Cock, PJA et al. The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants. Nuc. Acids Res. 2010 38(6):1767-1771). The nice thing with the conversion is that the sanger format does not list the read info for both the read sequence and read quality values so the files are significantly smaller.

Step 1 - Download the following source files and patches

UPDATE - This step is no longer necessary as recent versions of BWA allow you to do this on the fly during alignment using the -I option. I'll leave this here in case people need a way to process illumina 1.3-1.7 fastq into sanger format. Watch the version of casava used by your core in the coming year as the new illumina 1.8 pipeline will output files in sanger format so conversion will no longer be needed.

A) Maq (http://sourceforge.net/projects/maq/) *TO REALLY CONFUSE YOU THE MAIN DOWNLOAD IS ACTUALLY BWA, argh...)

NOTE: I only downloaded this to use the ill2sanger command to convert Illumina 1.3+ fastq files (ie. s_1_sequence.txt) to Sanger fastq format. Other options exist BioPerl, BioPython but I couldn't figure them out

- click on "View all files"
- click on "maq" folder
- click on newest version "0.7.1" (ASSUMPTION: This should be a long standing version since it appears that the development of Maq is dead with Heng Li releasing BWA)
- click on the download file "maq-0.7.1.tar.bz2"

Now download the ill2sanger patch that is needed to convert illumina 1.3+ fastq files to Sanger fastq

- click on "Develop" tab
- click on "Tracker" tab and select the "Patches" dropdown menu

NOTE: There are two versions to download the historic one by "daweonline" and a new alternative by "joelmartin" I used the original patch as I could find install instructions on Seqanswers

- click on ID 2841164 "illumina to sanger conversion"
- click download to download "maq-ill2sanger.patch" (current version submitted 2009-08-20)

Now Compile, patch, and re-compile the application

- move "maq-0.7.1.tar.bz2" to your "NGS/ApplicationDownloads" folder
- move "maq-ill2sanger.patch" to your "NGS" folder
- double click "maq-0.7.1.tar.bz2" to decompress the file
- open up "Terminal"
- navigate to the decompressed folder cd Documents/NGS/ApplicationDownloads/maq-0.7.1
- compile as per option 2 in Maq Manual (Release 0.5.0) (http://maq.sourceforge.net/maq-man.shtml)
> enter the following command "make -f Makefile.generic"
- a bunch of "Stuff" will come up, check to ensure no errors are listed!
- apply the patch as follows:
> step back one directory "cd .."
> run "ls" command to ensure the current directory contains folder "maq-0.7.1" and the patch "maq-ill2sanger.patch"
> install the patch with the following command "cd maq-0.7.1; patch -p1 < ../maq-ill2sanger.patch" found on seqanswers (http://seqanswers.com/forums/showthread.php?t=2499)
* You should get the following messages : patching file fastq2bfq.c
patching file main.c
patching file main.h
- recompile the maq application
> "make -f Makefile.generic"

Now check to see if the conversion patch was successful
> enter "./maq"

***This should bring up a window with the maq command options, check that ill2sanger is available under the "Format Converting" section***

B) BWA (http://sourceforge.net/projects/bio-bwa/files/)

- click on the download link for the newest version "bwa-0.5.9"
- download the file "bwa-0.5.9.tar.bz2"
- move "bwa-0.5.9.tar.bz2" to your "NGS/ApplicationDownloads" folder
- double click "bwa-0.5.9.tar.bz2" to decompress the file
- open up "Terminal"
- navigate to the decompressed folder "cd Documents/NGS/ApplicationDownloads/bwa-0.5.9"
- compile with make command
> enter the following command "make" (Really its that simple, this command line stuffs not that scary)
- a bunch of "Stuff" will come up, check to ensure no errors are listed!

Now check to see if the install was successful
> enter "./bwa"

***This should bring up a window with the bwa command options***

C) SAMtools (http://sourceforge.net/projects/samtools/)

- click on the download link for the newest version "samtools-0.1.12a"
- download the file "samtools-0.1.12a.tar.bz2"
- move "samtools-0.1.12a.tar.bz2" to your "NGS/ApplicationDownloads" folder
- double click "samtools-0.1.12a.tar.bz2" to decompress the file
- open up "Terminal"
- navigate to the decompressed folder "cd Documents/NGS/ApplicationDownloads/samtools-0.1.12a"
- compile with make command
> enter the following command "make" (Really its that simple)
- a bunch of "Stuff" will come up, check to ensure no errors are listed!

Now check to see if the install was successful
> enter "./samtools"

***This should bring up a window with the samtools command options***

The next step is to make each application executable as per the previous post options and we are just about ready to go.

Since we will use BWA we need to download the reference genomes to align against. The simplest place to get the data seems to be ensembl (http://www.ensembl.org/info/data/ftp/index.html) but we run into a problem with the full human genome file (Homo_sapiens.GRCh37.57.dna.toplevel.fa.gz) as it exceeds the maximum character length allowed by the bwa index command. To get around this problem if you want to use a GRCh37/hg19 genome version the best option seems to be the 1000 genomes version (ftp://ftp.sanger.ac.uk/pub/1000genom...ect_reference/) file (human_g1k_v37.fasta.gz). Copy the human_g1k_v37.fasta.gz file to the NGS/RefGenomes folder and then decompress it by double clicking on the file.

Test Question: Do you know the difference between UCSC mapping versus NCBI/Ensembl....Hint: Its a difference of 0 and 1 but it can really ruin your day when the commercial software manufacture doesn't know the difference!!

Next step, making the computer do some of the work

See next post,

Jonathan

[KevinLam]

I don't understand why you say there are problems with colorspace? Aligners (like BFAST) will convert your csfasta/qual files into FASTQ, will align your data sensitively, output to the SAM format, and then any SNP caller can be used without modification. Beyond specifying one command line option (to say the data is color space) during alignment there is no difference between Illumina/SOLiD (basespace/colorspace) data in terms of processing. It's the same workflow. Also the theoretical and practical benefits of colorspace (low false discovery rate) are rarely mentioned.

Sorry, just a slight pet peeve from a very happy SOLiD user.

Nils

Hi Nils,
No offence meant! It's all in good fun.
by problems I think I meant it more as caveats that you should watch for.

firstly it seems terribly important to understand dual base encoding but actually you just need an overview.

2ndly you are stuck with color space aware progs unless you wanna throw the benefits of colorspace away by direct conversion to base space and risk 3' ends being wrongly converted.
for de novo assembly with velvet you have to double encode your file into a format that looks exactly like 25 bp basespace fasta files. which can be misleading if someone else comes across the file and doesn't read the documentation you left there.

and it doesn't help that ABI's documentation for their software rarely exceeds 3 pages in pdf.

Other than that I am nearly a happy SOLiD user as you

[nilshomer]

rofl I shld do a version with SOLiD data with the rainbow assorted myriad of problems with colorspace. Good Post!

I don't understand why you say there are problems with colorspace? Aligners (like BFAST) will convert your csfasta/qual files into FASTQ, will align your data sensitively, output to the SAM format, and then any SNP caller can be used without modification. Beyond specifying one command line option (to say the data is color space) during alignment there is no difference between Illumina/SOLiD (basespace/colorspace) data in terms of processing. It's the same workflow. Also the theoretical and practical benefits of colorspace (low false discovery rate) are rarely mentioned.

Sorry, just a slight pet peeve from a very happy SOLiD user.

Nils

[KevinLam]

rofl I shld do a version with SOLiD data with the rainbow assorted myriad of problems with colorspace. Good Post!