Text Manipulation via Bash commands and Blast

Text Manipulation via Bash commands and Blast

A very common problem in bioinformatics is how to efficiently and effectively pull the information you want from large, sprawling data files. You could search online for pre-made Perl scripts or write your own Perl script to do this, but often it can be time-consuming and doesn't have a good cost:benefit ratio. Alternatively, you can learn how to do text manipulation using the Terminal, with sed commands and using MPIBlast. For instance, you can do text manipulation, using a blast database to pull your data. The idea is that you format a database, that contains the data in full. Then you blast, or send a query against the database. The query will have a list of information you want mpiblast to pull from the full dataset and output to a file. The query file is a text file that contains information identifying what sequences, etc. you want to pull out.

MPIBlast is an important and useful program to have. It is also a real pain to install on your computer.
To see if its installed already, open the terminal and type which mpiblast. I did this and got back:

/usr/local/bin/mpiblast

This tells me I have it installed somewhere on my Mac. This location just points to the program, but you want to see where the folder is, do a search. If you are lucky and have MPIBlast installed, then make sure that MPIBlast is in your path, but adding it to your bash profile. Your bash profile is usually in your home folder. If you are using Ubuntu it is called .bash_src. You can locate this by opening the terminal and cd into your home folder. Type ls -a to see all the hidden files, such as ones beginning with a '.' On a Mac:

cd /Users/username
ls -a
emacs .bash_profile
Add to the end:
PATH=$PATH\:/Applications/mpiblast-1.6.0/ncbi/bin ; export PATH
control+x+c
y     (to save)
source ~/.bash_profile  (reloads your bash profile while your Terminal is still open).

This strategy of text manipulation is a multi-step process, which will result in a series of output files that progressively move towards the final output you want. I am going to take a sample task in which I need to pull certain data from a large text file. Thankfully, a colleague was here to show me how to do this in a quick and efficient manner. Otherwise, I could easily waste many hours trying to Google ways to text manipulate. The sed command can be used in the Terminal to quickly and easily manipulate text. (There are tutorials online on how to use sed commands). You can either print to the screen or use the > fileoutput.txt to output the data to a file.

I have two major steps to do: reformat both the input file (containing the data I want to pull) and the fasta file. Then I can build a database. First off, I have a large file containing repeat elements, named consensi.fa.classified. I want to pull the sequences for the repeat elements labeled "Unknown." I have attached this sample file to this post.

The entries in this file look like:

>rnd-2_family-448#DNA ( Recon Family Size = 52 Final Multiple Alignment Size = 44 )
GGTGCATCTGCGGCAGGCATCCTCGGAGGCTGTGAGGTCTGGAGTTGTGG
GANTTGAAGTCCAAAACACCCGGAGGGCCGAAGTTTGCCCGTGCCTGATT
GACACTGTAGAATTAATGCAGTTTGACACCACTTTAACTGCCATGGCTCA
ATGCTATGGAATCCTGGGAGTTGTAGTTTGGTGAGGCACCAGCACTCTTT
GGCAGAGAAGGCTAAAGACCTCGTAAAACTACAACTCCCACGATTCCATA
GCATTGAGCCATGGCAGTTAAAGTGGTGTCAAACCGCATTAATTCTACAG
TGTAG

>rnd-2_family-1258#Unknown ( Recon Family Size = 45 Final Multiple Alignment Size = 43 )
AGGGGTGCTCACACTTTTTCAGCATGCGAGCTACTTATAAAACAAC

>rnd-2_family-160#LINE/RTE-BovB ( Recon Family Size = 53 Final Multiple Alignment Size = 41 )
GATTGATAAAGAGATAGATCACAGATTAGCAAAGGCATANAGTGCATTCG
GAAGGCTTCACAAAAGAGTCTGGAGCAACAAGCACCTGAGGCAAAGCACC
AAAATCAGTGTGTACAGAGCTATCGTACTGTCTATTCTCCTCTATGGGTC
TGAAACATGGGTCACCTATCGCCAACACCTACGACTCCTCGAACGCTTTC
ATCAGCGCTGTCTTCGCACAATCCTAAATATACACTGGACCGACTATGTG
ACGAATGTTGCTGTCCTTGAGCAAGCAGGGATCACCAGCGTTGAGGCCAT

Each entry is between a > and a >. The first line has information which acts as an identifier. It is difficult to create a script to pull the information from > to >, for every first line containing the word "Unknown." So what we will do is to first pull ONLY the first line for every entry with the word "Unknown."

I used the sed command:
sed -n '/Unknown/p' /Volumes/scratch/results/consensi.fa.classified > temp1.txt

I did this and got the output file temp1.txt, which looks like this:

rnd-2_family-1258#Unknown ( Recon Family Size = 45 Final Multiple Alignment Size = 43 )
rnd-2_family-144#Unknown ( Recon Family Size = 43 Final Multiple Alignment Size = 41 )
rnd-2_family-295#Unknown ( Recon Family Size = 44 Final Multiple Alignment Size = 40 )

When I try to run fastacmd, I get a long list of errors:

fastacmd -p F -o temp2.txt -i /Volumes/scratch/results/temp1.txt -d /Volumes/scratch/TEs/results/consensi.fa.classified

[fastacmd] ERROR: Entry "(" not found
[fastacmd] ERROR: Entry "recon" not found
[fastacmd] ERROR: Entry "family" not found
[fastacmd] ERROR: Entry "size" not found
[fastacmd] ERROR: Entry "=" not found
[fastacmd] ERROR: Entry "45" not found

Now, in sed, sometimes ( ) can cause confusion, as do spaces, so we want to get rid of them. 

We can do this by trading a space for a tab.

The command is:
sed 's/ /    /' /Volumes/scratch/results/consensi.fa.classified > newdb

You have to type control+v+TAB to create the tab in sed.
This sed command will take a space and turn it into a tab.
Tab-delineated files can be easy to work with. 


Note: Shortcut keys vary from system to system. In EMACS, you have to look into key binding to find out what the shortcuts are, or to set them to your desired key combinations and/or with meta-keys.


After the sed command, the output file temp3.txt looks like this:

rnd-2_family-1258#Unknown    ( Recon Family Size = 45 Final Multiple Alignment Size = 43 )
rnd-2_family-144#Unknown    ( Recon Family Size = 43 Final Multiple Alignment Size = 41 )
rnd-2_family-295#Unknown    ( Recon Family Size = 44 Final Multiple Alignment Size = 40 )
rnd-2_family-368#Unknown    ( Recon Family Size = 46 Final Multiple Alignment Size = 40 )

I want to get rid of the extra text: '    ( Recon Family Size = 46 Final Multiple Alignment Size = 40 )', since I only need the Unknown information.

So, I typed the following command to cut the text:
cut -f1 /Volumes/scratch/TEs/results/temp3.txt > temp4.txt

Then I got an output file called temp4.txt which looks like this:

rnd-2_family-1258#Unknown
rnd-2_family-144#Unknown
rnd-2_family-295#Unknown
rnd-2_family-368#Unknown

This particular cut command will delete everything after a tab. So if there are spaces or strange formatting, I do a replace all with a tab. (ie. change spaces behind Unknown to a tab behind Unknown).

Then you will need to use two commands: formatdb and fastcmd. If you are in doubt about how to use these tools, you can always type the command then help. Help varies from program to program, but in general its either '-h' or '--h' or '-help' and sometimes '--help'. For instance, to find out how to use mpiblast's formatdb command, type:

formatdb --help

This gives me a long list of options that I need/can specify. Using this knowledge, I typed the following commands:

formatdb -t repeat2.txt  -p false -o T -i /Volumes/scratch/results/finaldb

Similarly, to figure out how to use mpiblast's fastacmd command, I typed:

fastacmd --help
fastacmd -p F -o 5Jul_Unknown_Elements.txt -i /Volumes/scratch/TEs/results/temp4.txt -d  /Volumes/scratch/results/finaldb


I get an output file that, for some reason, has two things added to it now that I don't want:

>lcl|rnd-2_family-1258#Unknown No definition line found
AGGGGTGCTCACACTTTTTCAGCATGCGAGCTACTTATAAAACAAC
>lcl|rnd-2_family-144#Unknown No definition line found
AATGATGTAGACCAGGCATGGGCAAACTTGGGCCCTCCAGGTGTTTTGGACTTCAGGCGTTTTGGACTTCAANTCCCAGA
ATTCCTAGCAGCCTGCAGGCTGTTAGGAATTGTGGGAGTTGAAGTCCAAAACACCTGGAGGGCCGAAGTTTGCCCATGCC
TGATCTACACTGTAGAATTAATGCAGTTTGACACCACTTTAACTGCCATGGCTCAACGCTATGGAATCCTGGGAGTTGTA
GTTTGGCGAGGCACCAGCACTCTTTGGCAGAGAAGGCTAAAGAAAACTACAACTCCCATGATTCCATAGCAGTTGAAGTG
NTGTCAAACTGCATCAATTGTACAGNTGTGT
>lcl|rnd-2_family-295#Unknown No definition line found
AGAACATAAGAACATAAGAAGAGCCATGCTGAATCGGGCCTCAGCCC
>lcl|rnd-2_family-368#Unknown No definition line found
GTACACTGTAGAATTAATGCAGTTTGACACCACTTTAACTGCCATGGCTCAATGCTATGGAATCATGGGAGTTGTAGTTT
TACAAGGTCTTCTCTGCCAAAGAGTGCTGGTGCCTCGCCAAACTACAANTCCCAGGATCCCATAGCATTGAGCCATGGCA
GTTAAAGTGGTGTCAAACTGCATTAATTCTACAGTGTAG


In the text editor, I do a simple find and replace for '>lcl|' to nothing and 'No definition line found' to nothing.
This gets rid of the extra text and gives me my final output:

rnd-2_family-1258#Unknown
AGGGGTGCTCACACTTTTTCAGCATGCGAGCTACTTATAAAACAAC
rnd-2_family-144#Unknown
AATGATGTAGACCAGGCATGGGCAAACTTGGGCCCTCCAGGTGTTTTGGACTTCAGGCGTTTTGGACTTCAANTCCCAGA
ATTCCTAGCAGCCTGCAGGCTGTTAGGAATTGTGGGAGTTGAAGTCCAAAACACCTGGAGGGCCGAAGTTTGCCCATGCC
TGATCTACACTGTAGAATTAATGCAGTTTGACACCACTTTAACTGCCATGGCTCAACGCTATGGAATCCTGGGAGTTGTA
GTTTGGCGAGGCACCAGCACTCTTTGGCAGAGAAGGCTAAAGAAAACTACAACTCCCATGATTCCATAGCAGTTGAAGTG
NTGTCAAACTGCATCAATTGTACAGNTGTGT
rnd-2_family-295#Unknown
AGAACATAAGAACATAAGAAGAGCCATGCTGAATCGGGCCTCAGCCC
rnd-2_family-368#Unknown
GTACACTGTAGAATTAATGCAGTTTGACACCACTTTAACTGCCATGGCTCAATGCTATGGAATCATGGGAGTTGTAGTTT
TACAAGGTCTTCTCTGCCAAAGAGTGCTGGTGCCTCGCCAAACTACAANTCCCAGGATCCCATAGCATTGAGCCATGGCA
GTTAAAGTGGTGTCAAACTGCATTAATTCTACAGTGTAG
rnd-2_family-267#Unknown
CTTGGCNTTNCTTCGCGAACGAAGATTCTTAGGAAGGACTCATCCCACGCTTTCTACAAGCGCGTTGATGATTCGGAGGT
TGGATGTGTCCTTCGAGCCTGTGCAATGGATTTTTCTGGTGGAGCGCAGCGTGCACAGAACTGGCCTCACCCTTTAAACC
AAAGGTTCGTCTGCCGGGGCCTAGCAAGCTTGGACGGTGGCAGCGAGGTCCTCAGGTCGTAGGTTTGGTTCGAGT

Another great tool is cat! For instance, I had two large files about ~1 GB each that I needed to be in one file. With smaller files, I could copy and paste the text of one file into the text of the other, etc. But this is not really possible with files this large. Using the bash command cat (concatenate) is the easy answer. 

cat    file1.txt    file2.txt    > /Users/username/Desktop/appended_file.txt

Now, I have an appended file with all the data I need.

I am just learning how to do work with these tools, but I hope this helps you out!

Saguaro HPC Computing

Supercomputing at ASU:

Using Saguaro

You will need to email the staff at AC2C Computing center to have a Saguro account added. Once you have your account you can log on via SSH and a Terminal. Mac and Linux-based systems (like Ubuntu) have a Terminal built in. Windows is more complicated, but you can download some programs for ssh. Also, you can download VMWare and run a virtual drive with Ubuntu from a Windows operating system. This takes a bit of time, but I recommend it, because its nice to use. 

To log onto Saguaro, open up a Terminal and type:

ssh username@saguaro.fulton.asu.edu

password ************

To find your way around you can type cd and ls and cd .. to navigate and peek at the files in the folder directory. You will probably want to set up new folders to begin with. To make a new folder type: mkdir new_folder.

I recommend using lower case letters and avoid using spaces, instead_using_underscores. To find hidden files, like your .bash_profile file on your home Saguro folder, type:  ls -a.

To set permissions on folders or files:     

ls -l      chmod 700 file    

You will see something like:

rwxr-xr-x username user bio.pl

An example of the permissions encoded might be rwxr-xr-x. What this means is that the owner can do anything with the file, but group owners and the rest of the world can only read or execute it. (w-write,xr-execute and read,x-execute).

There is an easy way to do this and a harder way. Just memorize the permissions codes. Read access is the number 4, write is 2, and execute is 1 (4+2+1=7). So for you and you only, 7 is the permission you want for yourself.

chmod 700   bio.txt      Only you can read, write to, or execute bio.txt

chmod 777   bio.txt      Everybody can read, write to, or execute bio.txt

chmod 744   bio.txt      Only you can read, write to, or execute bio.txt. Everybody can read bio.txt

If you want to set permissions for an entire folder type:

chmod  -R   700   /home/username/folder_name OR if you are already in the directory where the folder is

chmod  -R   700   folder_name

The hard way is to learn the permissions coding system, which I will go ahead and tell you: Every permissions mode has a corresponding code number, and one number that corresponds to any mode. Every one of the three digits on the mode number corresponds to one of the three permission triplets. (u, g and o) Every permission bit in a triplet corresponds to a value: 4 for r, 2 for w, 1 for x. If the permission bit you add this value to the number of the permission triplet. So if a file has rwxr-xr-x permissions we do the following calculation:

Triplet for u: rwx => 4 + 2 + 1 = 7  (user)

Triplet for g: r-x => 4 + 0 + 1 = 5  (group owners)

Triplet for o: r-x => 4 + 0 + 1 = 5  (world)

That gives: 755. In other words, 755 is a code way to say 'Other people can read or run this file, but only I should be able to modify it' and 777 means 'everyone has full access to this file.'

#	Permission	rwx
7	full	111
6	read and write	110
5	read and execute	101
4	read only	100
3	write and execute	011
2	write only	010
1	execute only	001
0	none	000

It's a good idea to get used setting permissions correctly, so I recommend making a new directory, and a few new files to get used to it.

mkdir test

cd test

emacs test1.txt

ctrl+x and ctrl+c

cp test1.txt test2.txt

cp test2.txt text3.txt

ls -l

chmod 700 test1.txt

chmod 755 test2.txt

ls -l

cd ..

chmod  -R   700   test

cd test

ls -l

On Saguro, you want to exclude other users from accessing your files openly. Become familiar with the read-write-execute permissions. But BE CAREFUL!!!!! because if you set the permission wrong, you may lock yourself out of your own files!!!

Moving Large Files:

When you are transferring large files, you will need to use Datamover. If you do not, you will be warned by the computing staff about locking up the logon nodes. You will need to go to system preferences and find out your IP and your logon name,

To download stuff from Saguaro back to your local computer:

First, ssh into Saguaro. Then type:

scp  name@datamover1.fulton.asu.edu:/home/username/file.txt  username@111.111.11:/Users/username/Desktop/

name@datamover1.fulton.asu.edu's password:

Password: **********

file.txt                                                               100%  385KB 385.0KB/s   00:00

Use -r if you want to move a whole folder. Alternatively, you can compress or tar the file before moving it.

scp -r   name@datamover1.fulton.asu.edu:/scratch/username   /Volumes/scratch/

To download stuff from the local computer up to Saguaro:

From your local Terminal

scp /Users/username/Desktop/file.txt   username@datamover1.fulton.asu.edu:/home/username

Password: ******

file.txt                                                               100%  385KB 385.0KB/s   00:00

Don't forget that instead of typing out a long file path name, you can drag the file over to Terminal and it will record the path, after your cursor location. This only works locally, of course.

Okay, there are two spaces where you will store your stuff on Saguaro. Use a : after datamover to specify where you want the files to go.

1) At your home folder

    ex. /home/username/progs

2) On your scratch space

    ex. /scratch/username/25Jun

You don't want to datamove huge files to your home drive, since space is limited. You can only store large drives on your scratch space for a few months at a time. You will usually receive a notice every 30 days that the scratch space will be purged. 

To submit jobs to saguaro, you will need a script. This is just a text file that will tell the job how to run. You can literally copy this into a text file and save and use on Saguaro.   

script.txt:

--------------------------------------------------------------

#!/bin/bash

#PBS -l nodes=12

#PBS -l walltime=96:00:00

#PBS -l nodeset=ANYOF:FEATURE:westmere

#PBS -j oe

#PBS -o /scratch/username/aapl49

module load python

cd /scratch/username2/alignments/Acar_genome

bowtie2-build Anolis_genome.fa Anolis_genome

exit

----------------------------------------------------------------

The -o is the output directory. The nodeset can be set to other machines, or left to random. Check the Saguro .pdf for more info. Sometimes you will need modules (preinstalled software on Saguaro). The max time for a job is 96 hours. Make sure that you cd into the correct folder, or if your files aren't in the same folder that you specify the full path name to each file.

Often you will need to edit these scripts right on Saguro. Do this with the emacs command. Enter the directory of the file you want to edit. Then type:

emacs script.txt

Emacs is important to use, so its best to Google Emacs commands and get used to it! To save a quit type control-x and control-c at the same time. To suspend the shell, you can type control+z. To save as you can control+w, then typing the name of the file after the file path before hitting enter. Also, the terminal will only refresh if all windows are closed, then you open it again. If you want to refresh it with your window still open, type: source ~/.bash_profile.

When you are ready to submit a job, type:

qsub script.txt

To check on your job you can type either:

qstat -a

watch qstat -a

You can note the job ID number from qstat -a.

To watch an activity monitor ranked by CPU used:

top -o cpu

Often, your job will be in the queue and may take a while. To see what jobs are running before you, you can type:

showq 

pbsnodes  (also gives you information about what nodes are being used)

Often you will need to check if a job is running correctly. You can collect information from:

qpeek jobnumber

qmem job number

For example:

qmem 5102197

+---------+---------------------+---------------------+---------------------+

| Host    |   RAM (Used/Total)           Swap (Used/Total)   | Load Average |

+---------+---------------------+---------------------+---------------------+

|   s6-6  |   15.57GB / 15.67GB |       7.23GB / 7.81GB |    2.56, 1.76, 1.56 |

|  s21-10 |  751.73MB / 15.67GB |    30.93MB / 7.81GB |  0.00, 0.00, 0.00 |

+---------+---------------------+---------------------+---------------------+

If you see gigabytes of Swap being used, this is not a good sign and you should cancel the job and try to run it with more processors, nodes or on the very large SMP machine.

To run a job on the SMP (32-core) machine, in the script just specify:

#PBS   -q smp

To delete a job:

qdel jobnumber

ex. qdel 477551

To hold a job:

qhold jobnumber

To release a hold:

qrls jobnumber

SNPs with SamTools

SNPs with SamTools

These are kind of my messy notes on SNP bioinformatics. Its best to consult the manuals and/or documentation first. SeqAnswers also is a great source to find a Q&A forum of other bioinformaticians doing similar work. A screenshot is attached to show to expected file output.

Bowtie and Tophat are run on a shared resource HPC computer. You can install them on your home folder or you can use module load when you write your script. To see what programs are already on Saguaro, you can type module avail. On your local machine, you can run http://samtools.sourceforge.net/cns0.shtml SamTools and BCFTools (part of Samtools). You should also probably have the IGV Genome browser, as well. SamTools can be used to find variations. It is a suite for storing, manipulating and analyzing alignments, like Bowtie outputs Samtools has two formats- user friendly SAM and binary BAM output. Sam tools can be used to find SNPs in a Bowtie output file. Run bowtie to align the reads, -S for sam output.


Here is the Saguaro script: 

------------------------------------------------------------------
#!/bin/bash

#PBS -l nodes=12
#PBS -l walltime=96:00:00
#PBS -j oe
#PBS -o /scratch/username/aapl49

cd /scratch/username/
module load bowtie
tophat -r 139 -o aapl49 -p 12 /scratch/username/kmer56_q20 Sample_DNA.fastq  Sample_DNA.fastq

exit
-------------------------------------------------------------------


I got all the files back from Tophat and downloaded them to my local drive from Saguaro.

Don't forget to add SamTools to your path.
cd /home/username 
ls -a
open .bash_profile
PATH=$PATH\:/Users/username/progs; export PATH
PATH=$PATH\:/Users/username/progs/samtools118/bcftools ; export PATH

We want accepted_hits.BAM file
(.bed will be a tab-delineated file, with insertions and deletions)
Put Sam Tools in your path, otherwise go in folder and run ./samtools

General Pattern:
samtools index accepted_hits.bam
samtools view -b accepted_hits.bam chr2 > accepted_hits.chr2.bam
samtools sort accepted_hits.chr2.bam file.bam

Convert SAM into BAM to sort
samtools view -bS -o ec_snp.bam ec_snp.sam
Then the BAM file is sorted (preparing for SNP calling)
samtools sort ec_snp.bam ec_snp.sorted

Sort the BAM files
samtools sort [-no] [-m maxMem] <in.bam> <out.prefix>
./samtools sort /Volumes/scratch/aapl49/accepted_hits49.bam /Volumes/scratch/aapl49/accepted_hits49_sorted
./samtools sort /Volumes/scratch/aapl75/accepted_hits75.bam /Volumes/scratch/aapl75/accepted_hits75_sorted

My colleague says that I don't need to sort or index the files, as they come out of Tophat already sorted. 

What I do need is the .fasta file to be indexed into a .fai file


SamTools SNP calling procedure:
From a sorted BAM alignment, raw SNP and indel calls are acquired by:
samtools mpileup /Volumes/scratch/aapl49/accepted_hits49_sorted.bam > raw.49.pileup
samtools  mpileup /Volumes/scratch/aapl75/accepted_hits75_sorted.bam > raw.75.pileup

To extract a certain sequence, faidx will index the file and create <ref.fasta>.fai on the disk
samtools faidx /Volumes/scratch/kmer56_q20.fa [region1 [...]]
you can rename .fa files to .fasta (its just the same)

samtools faidx ref.fasta

If you get this error message: [fai_build_core] different line length in sequence 'scaffold50'. Segmentation fault. Use:
sed '/^$/d' myFile > tt
to reformat the linebreaks.

Learn how to reformat text files in the BASH using the SED command! (cut/sort/replace/find)
samtools faidx /Volumes/scratch/kmer56_q20.fasta
The > denotes outputting the file (instead of displaying it in bash).
On the MacOSX, the general text editor will default to saving files in a .rtf format. This is a pain. You can change the setting so the default is a .txt file. I highly recommend doing this to prevent create symbols/problems in running things later on.


My colleague Dr. Walter Eckalbar emailed me the following information:

samtools mpileup -uf <genome> <bam1> <bam2> | bcftools view -bvcg - > <var.raw.bcf>

#To make vcf file:
&& bcftools view <var.raw.bcf> | vcfutils.pl varFilter -D10 > <var.flt.vcf>

#To make index file
bgzip var.flt.vcf
tabix -p vcf var.flt.vcf.gz


samtools mpileup -uf <genome> <bam1> <bam2> | bcftools view -bvcg - > <var.raw.bcf>
samtools mpileup -uf /Volumes/scratch/kmer56_q20.fasta  /Volumes/scratch/aapl49/49snps.in.bam /Volumes/scratch/aaplgg70/gg70snps.in.bam /Volumes/scratch/aapl75/75snps.in.bam /Volumes/scratch/aapl20a/20asnps.in.bam  | bcftools view -bvcg - > /Volumes/scratch/49_75_gg70_20a_snps.var.raw.bcf
This has all 4 individuals for comparison.

#To make vcf file:
&& bcftools view <var.raw.bcf> | vcfutils.pl varFilter -D10 > <var.flt.vcf>
&& bcftools view kmer56_q20.49to75.var.raw.bcf | vcfutils.pl varFilter -D10 > kmer56_q20.49to75.var.flt.vcf
&& bcftools view /Volumes/scratch/kmer56_q20.49to75.var.raw.bcf  | vcfutils.pl varFilter -D10 > /Volumes/scratch/kmer56_q20.49to75.var.flt.vcf

#To make index file
bgzip var.flt.vcf
bgzip /Volumes/scratch/kmer56_q20.49to75.var.flt.vcf
tabix -p vcf var.flt.vcf.gz
tabix -p vcf /Volumes/scratch/kmer56_q20.49to75.var.flt.vcf.gz
then I get a kmer56_q20.49to75.var.flt.vcf.gz.tbi file


You will want the program IGV to browse the files then. This is a Java program made by the Broad Institute. 


Additional Notes:
Quality index is logarithmic:
Q10 means that there is a 10% error rate
Q20 0.01 1% error cutoff
Q30 is 0.001 chance of error

Fasta is a reference file.
Raw data to align are fastq files.


Running Bowtie will give .ebwt, this is the indexed reference file.
Then use .ebwt as an input for alignment in Topcoat.
Create new directories for each run in topcoat, because it will generate non-unique file names.

Tophat better for picking up indels between the species.

Reading the output file: Where are the SNPS? 

If this is one of the lines in your pileup file:
chr2 218208 G 30 ....................,......... +(+++)++++!#+++!+)+++!*****(#(

The 5th column with the dots and commas is the one of interest. Those are the "reads" that map to that position.
What you are looking for is letters instead of dots and commas, those are "mismatches" to the reference and in theory your SNPs (of course, you have to have plenty of mismatches to the reference for it to be considered a variant).

First, sort the BAM files

To Index or Not to Index?
samtools index /Volumes/scratch/aapl49/accepted_hits49_sorted.bam
samtools index /Volumes/scratch/aapl75/accepted_hits75_sorted.bam 

Index file <aln.bam>.bai will be created.

From this I got the file outputs:
/Volumes/scratch/aapl49/accepted_hits49_sorted.indexed.bam.bai
/Volumes/scratch/aapl75/accepted_hits75_sorted.indexed.bam.bai

To check help: samtools mpileup
samtools view /Volumes/scratch/aapl49/accepted_hits49_sorted.bam chr2:20,100,000-20,200,000
samtools view /Volumes/scratch/aapl75/accepted_hits75_sorted.bam chr2:20,100,000-20,200,000
samtools tview accepted_hits49_sorted.bam /Volumes/scratch/kmer56_q20.fa

samtools mpileup -C50 -gf ref.fasta -r chr3:1,000-2,000 in1.bam in2.bam
samtools tview aln.sorted.bam ref.fasta

It might be better just to call for a certain chromosome at a time. 


Retrieve and print stats in the index file. The output is a TAB delimited file with each line consisting of reference sequence name, sequence length, # mapped reads and # unmapped reads.
samtools idxstats /Volumes/scratch/aapl49/accepted_hits49_sorted.bam
samtools idxstats /Volumes/scratch/aapl75/accepted_hits75_sorted.bam

Even though i specified an output directory, the raw.pileup and raw49.pileup files were placed in the samtools directory

The resultant output should be further filtered by:
samtools.pl varFilter raw.pileup | awk '$6>=20' > final.pileup
perl /Users/username/progs/samtools118/samtools.pl varFilter raw.pileup | awk '$6>=20' > final.pileup
to rule out error-prone variant calls caused by factors not considered in the statistical model.