cd) to the linuxTutorial folder, and do all the processing inside:
[uzi@quince-srv2 ~/]$ cp -r /home/opt/MScBioinformatics/linuxTutorial . [uzi@quince-srv2 ~/]$ cd linuxTutorial [uzi@quince-srv2 ~/linuxTutorial]$I have deliberately chosen
Awk in the exercises as it is a language in itself and is used more often to manipulate NGS data as compared to the other command line tools such as grep, sed, perl etc. Furthermore, having a command on awk will make it easier to understand advanced tutorials such as Illumina Amplicons Processing Workflow.
Linux, we use a shell that is a program that takes your commands from the keyboard and gives them to the operating system. Most Linux systems utilize Bourne Again SHell (bash), but there are several additional shell programs on a typical Linux system such as ksh, tcsh, and zsh. To see which shell you are using, type
[uzi@quince-srv2 ~/linuxTutorial]$ echo $SHELL /bin/bashTo see where you are in the file system:
[uzi@quince-srv2 ~/linuxTutorial]$ pwd /home/uzi/linuxTutorialList the files in the current directory:
[uzi@quince-srv2 ~/linuxTutorial]$ ls dataNow try different commands from the sheet given below:
ls — list items in current directory
ls -l — list items in current directory and show in long format to see perimissions, size, and modification date
ls -a — list all items in current directory, including hidden files
ls -F — list all items in current directory and show directories with a slash and executables with a star
ls dir — list all items in directory dir
cd dir — change directory to dir
cd .. — go up one directory
cd / — go to the root directory
cd ~ — go to to your home directory
cd - — go to the last directory you were just in
pwd — show present working directory
mkdir dir — make directory dir
rm file — remove file
rm -r dir — remove directory dir recursively
cp file1 file2 — copy file1 to file2
cp -r dir1 dir2 — copy directory dir1 to dir2 recursively
mv file1 file2 — move (rename) file1 to file2
ln -s file link — create symbolic link to file
touch file — create or update file
cat file — output the contents of file
less file — view file with page navigation
head file — output the first 10 lines of file
tail file — output the last 10 lines of file
tail -f file — output the contents of file as it grows, starting with the last 10 lines
vim file — edit file
alias name 'command' — create an alias for a command
shutdown — shut down machine
reboot — restart machine
date — show the current date and time
whoami — who you are logged in as
finger user — display information about user
man command — show the manual for command
df — show disk usage
du — show directory space usage
free — show memory and swap usage
whereis app — show possible locations of app
which app — show which app will be run by default
ps — display your currently active processes
top — display all running processes
kill pid — kill process id pid
kill -9 pid — force kill process id pid
ls -l — list items in current directory and show permissions
chmod ugo file — change permissions of file to ugo - u is the user's permissions, g is the group's permissions, and o is everyone else's permissions. The values of u, g, and o can be any number between 0 and 7.
7 — full permissions
6 — read and write only
5 — read and execute only
4 — read only
3 — write and execute only
2 — write only
1 — execute only
0 — no permissions
chmod 600 file — you can read and write - good for files
chmod 700 file — you can read, write, and execute - good for scripts
chmod 644 file — you can read and write, and everyone else can only read - good for web pages
chmod 755 file — you can read, write, and execute, and everyone else can read and execute - good for programs that you want to share
wget file — download a file
curl file — download a file
scp user@host:file dir — secure copy a file from remote server to the dir directory on your machine
scp file user@host:dir — secure copy a file from your machine to the dir directory on a remote server
scp -r user@host:dir dir — secure copy the directory dir from remote server to the directory dir on your machine
ssh user@host — connect to host as user
ssh -p port user@host — connect to host on port as user
ssh-copy-id user@host — add your key to host for user to enable a keyed or passwordless login
ping host — ping host and output results
whois domain — get information for domain
dig domain — get DNS information for domain
dig -x host — reverse lookup host
lsof -i tcp:1337 — list all processes running on port 1337
grep pattern files — search for pattern in files
grep -r pattern dir — search recursively for pattern in dir
grep -rn pattern dir — search recursively for pattern in dir and show the line number found
grep -r pattern dir --include='*.ext — search recursively for pattern in dir and only search in files with .ext extension
command | grep pattern — search for pattern in the output of command
find file — find all instances of file in real system
locate file — find all instances of file using indexed database built from the updatedb command. Much faster than find
sed -i 's/day/night/g' file — find all occurrences of day in a file and replace them with night - s means substitude and g means global - sed also supports regular expressions
tar cf file.tar files — create a tar named file.tar containing files
tar xf file.tar — extract the files from file.tar
tar czf file.tar.gz files — create a tar with Gzip compression
tar xzf file.tar.gz — extract a tar using Gzip
gzip file — compresses file and renames it to file.gz
gzip -d file.gz — decompresses file.gz back to file
ctrl+a — move cursor to beginning of line
ctrl+f — move cursor to end of line
alt+f — move cursor forward 1 word
alt+b — move cursor backward 1 word
http://cheatsheetworld.com/programming/unix-linux-cheat-sheet/
Awk is a programming language which allows easy manipulation of structured data and is mostly used for pattern scanning and processing. It searches one or more files to see if they contain lines that match with the specified patterns and then perform associated actions.
The basic syntax is:
awk '/pattern1/ {Actions}
/pattern2/ {Actions}' file
The working of Awk is as follows
Awk reads the input files one line at a time.Awk performs the given actions for each line of the input.data.tsv with the following contents:
$ cat data/test.tsv blah_C1 ACTGTCTGTCACTGTGTTGTGATGTTGTGTGTG blah_C2 ACTTTATATATT blah_C3 ACTTATATATATATA blah_C4 ACTTATATATATATA blah_C5 ACTTTATATATTBy default
Awk prints every line from the file.
$ awk '{print;}' data/test.tsv
blah_C1 ACTGTCTGTCACTGTGTTGTGATGTTGTGTGTG
blah_C2 ACTTTATATATT
blah_C3 ACTTATATATATATA
blah_C4 ACTTATATATATATA
blah_C5 ACTTTATATATT
We print the line which matches the pattern blah_C3
$ awk '/blah_C3/' data/test.tsv blah_C3 ACTTATATATATATA
Awk has number of builtin variables. For each record i.e line, it splits the record delimited by whitespace character by default and stores it in the $n variables.
If the line has 5 words, it will be stored in $1, $2, $3, $4 and $5. $0 represents the whole line. NF is a builtin variable which represents the total number of fields in a record.
$ awk '{print $1","$2;}' data/test.tsv
blah_C1,ACTGTCTGTCACTGTGTTGTGATGTTGTGTGTG
blah_C2,ACTTTATATATT
blah_C3,ACTTATATATATATA
blah_C4,ACTTATATATATATA
blah_C5,ACTTTATATATT
$ awk '{print $1","$NF;}' data/test.tsv
blah_C1,ACTGTCTGTCACTGTGTTGTGATGTTGTGTGTG
blah_C2,ACTTTATATATT
blah_C3,ACTTATATATATATA
blah_C4,ACTTATATATATATA
blah_C5,ACTTTATATATT
Awk has two important patterns which are specified by the keyword called BEGIN and END. The syntax is as follows:
BEGIN { Actions before reading the file}
{Actions for everyline in the file}
END { Actions after reading the file }
For example,
$ awk 'BEGIN{print "Header,Sequence"}{print $1","$2;}END{print "-------"}' data/test.tsv
Header,Sequence
blah_C1,ACTGTCTGTCACTGTGTTGTGATGTTGTGTGTG
blah_C2,ACTTTATATATT
blah_C3,ACTTATATATATATA
blah_C4,ACTTATATATATATA
blah_C5,ACTTTATATATT
-------
We can also use the concept of a conditional operator in print statement of the form print CONDITION ? PRINT_IF_TRUE_TEXT : PRINT_IF_FALSE_TEXT. For example, in the code below, we identify sequences with lengths > 14:
$ awk '{print (length($2)>14) ? $0">14" : $0"<=14";}' data/test.tsv
blah_C1 ACTGTCTGTCACTGTGTTGTGATGTTGTGTGTG>14
blah_C2 ACTTTATATATT<=14
blah_C3 ACTTATATATATATA>14
blah_C4 ACTTATATATATATA>14
blah_C5 ACTTTATATATT<=14
We can also use 1 after the last block {} to print everything (1 is a shorthand notation for {print $0} which becomes {print} as without any argument print will print $0 by default), and within this block, we can change $0, for example to assign the first field to $0 for third line (NR==3), we can use:
$ awk 'NR==3{$0=$1}1' data/test.tsv
blah_C1 ACTGTCTGTCACTGTGTTGTGATGTTGTGTGTG
blah_C2 ACTTTATATATT
blah_C3
blah_C4 ACTTATATATATATA
blah_C5 ACTTTATATATT
You can have as many blocks as you want and they will be executed on each line in the order they appear, for example, if we want to print $1 three times (here we are using printf instead of print as the former doesn't put end-of-line character),
$ awk '{printf $1"\t"}{printf $1"\t"}{print $1}' data/test.tsv
blah_C1 blah_C1 blah_C1
blah_C2 blah_C2 blah_C2
blah_C3 blah_C3 blah_C3
blah_C4 blah_C4 blah_C4
blah_C5 blah_C5 blah_C5
Although, we can also skip executing later blocks for a given line by using next keyword:
$ awk '{printf $1"\t"}NR==3{print "";next}{print $1}' data/test.tsv
blah_C1 blah_C1
blah_C2 blah_C2
blah_C3
blah_C4 blah_C4
blah_C5 blah_C5
$ awk 'NR==3{print "";next}{printf $1"\t"}{print $1}' data/test.tsv
blah_C1 blah_C1
blah_C2 blah_C2
blah_C4 blah_C4
blah_C5 blah_C5
You can also use getline to load the contents of another file in addition to the one you are reading, for example, in the statement given below, the while loop will load each line from test.tsv into k until no more lines are to be read:
$ awk 'BEGIN{while((getline k <"data/test.tsv")>0) print "BEGIN:"k}{print}' data/test.tsv
BEGIN:blah_C1 ACTGTCTGTCACTGTGTTGTGATGTTGTGTGTG
BEGIN:blah_C2 ACTTTATATATT
BEGIN:blah_C3 ACTTATATATATATA
BEGIN:blah_C4 ACTTATATATATATA
BEGIN:blah_C5 ACTTTATATATT
blah_C1 ACTGTCTGTCACTGTGTTGTGATGTTGTGTGTG
blah_C2 ACTTTATATATT
blah_C3 ACTTATATATATATA
blah_C4 ACTTATATATATATA
blah_C5 ACTTTATATATT
You can also store data in the memory with the syntax VARIABLE_NAME[KEY]=VALUE which you can later use through for (INDEX in VARIABLE_NAME) command:
$ awk '{i[$1]=1}END{for (j in i) print j"<="i[j]}' data/test.tsv
blah_C1<=1
blah_C2<=1
blah_C3<=1
blah_C4<=1
blah_C5<=1
[uzi@quince-srv2 ~/linuxTutorial]$ cat data/test.fa >blah_C1 ACTGTCTGTC ACTGTGTTGTG ATGTTGTGTGTG >blah_C2 ACTTTATATATT >blah_C3 ACTTATATATATATA >blah_C4 ACTTATATATATATA >blah_C5 ACTTTATATATTand an IDs file:
[uzi@quince-srv2 ~/linuxTutorial]$ cat data/IDs.txt blah_C4 blah_C5After looking at the file, it is immediately clear that the sequences may span multiple lines (for example, for
blah_C1). If we want to match an ID, we can first linearize the file by using the conditional operator as discussed above to have the delimited information of each sequence in one line, and then make logic to perform further functionality on each line later.
Our logic is that for lines that contain header information /^>/ we can do something differently, and for other lines we use printf to remove new line character:
[uzi@quince-srv2 ~/linuxTutorial]$ awk '{printf /^>/ ? $0 : $0}' data/test.fa >blah_C1ACTGTCTGTCACTGTGTTGTGATGTTGTGTGTG>blah_C2ACTTTATATATT>blah_C3ACTTATATATATATA>blah_C4ACTTATATATATATA>blah_C5ACTTTATATATTWe can then put each sequence on a separate line and also put a tab character (
"\t") between the header and the sequence:
[uzi@quince-srv2 ~/linuxTutorial]$ awk '{printf /^>/ ? "\n"$0 : $0}' data/test.fa >blah_C1ACTGTCTGTCACTGTGTTGTGATGTTGTGTGTG >blah_C2ACTTTATATATT >blah_C3ACTTATATATATATA >blah_C4ACTTATATATATATA >blah_C5ACTTTATATATT[uzi@quince-srv2 ~/linuxTutorial]$ awk '{printf /^>/ ? "\n"$0"\t" : $0}' data/test.fa >blah_C1 ACTGTCTGTCACTGTGTTGTGATGTTGTGTGTG >blah_C2 ACTTTATATATT >blah_C3 ACTTATATATATATA >blah_C4 ACTTATATATATATA >blah_C5 ACTTTATATATTWe can then use
NR==1 block to stop printing a new line character before the first header (as you can see there is an empty space) and use next to ignore the later block:
[uzi@quince-srv2 ~/linuxTutorial]$ awk 'NR==1{printf $0"\t";next}{printf /^>/ ? "\n"$0"\t" : $0}' data/test.fa >blah_C1 ACTGTCTGTCACTGTGTTGTGATGTTGTGTGTG >blah_C2 ACTTTATATATT >blah_C3 ACTTATATATATATA >blah_C4 ACTTATATATATATA >blah_C5 ACTTTATATATTWe can then pipe this stream to another
awk statement using "\t" as a delimeter (which you can specify using -F) and use gsub to remove > from the start of each line since our IDs file doesn't contain that character:
[uzi@quince-srv2 ~/linuxTutorial]$ awk 'NR==1{printf $0"\t";next}{printf /^>/ ? "\n"$0"\t" : $0}' data/test.fa | awk -F"\t" '{gsub("^>","",$0);print $0}' blah_C1 ACTGTCTGTCACTGTGTTGTGATGTTGTGTGTG blah_C2 ACTTTATATATT blah_C3 ACTTATATATATATA blah_C4 ACTTATATATATATA blah_C5 ACTTTATATATTNow we load the
IDs.txt file in the BEGIN block, store the IDs in the memory, and in the stream if the first field ($1) matches the ID stored in the memory, we output the formatted record:
[uzi@quince-srv2 ~/linuxTutorial/data]$ awk 'NR==1{printf $0"\t";next}{printf /^>/ ? "\n"$0"\t" : $0}' data/test.fa | awk -F"\t" 'BEGIN{while((getline k < "data/IDs.txt")>0)i[k]=1}{gsub("^>","",$0); if(i[$1]){print ">"$1"\n"$2}}' >blah_C4 ACTTATATATATATA >blah_C5 ACTTTATATATTWith
Bioawk it is much simpler as you don't have to linearize the FASTA file as the record boundaries are the complete sequence boundaries and not lines:
[uzi@quince-srv2 ~/linuxTutorial/data]$ bioawk -cfastx 'BEGIN{while((getline k <"data/IDs.txt")>0)i[k]=1}{if(i[$name])print ">"$name"\n"$seq}' data/test.fa
>blah_C4
ACTTATATATATATA
>blah_C5
ACTTTATATATT
-c, with the field names as follows (you can use the column pairs alternatively):
bed: $1:$chrom $2:$start $3:$end $4:$name $5:$score $6:$strand $7:$thickstart $8:$thickend $9:$rgb $10:$blockcount $11:$blocksizes $12:$blockstarts sam: $1:$qname $2:$flag $3:$rname $4:$pos $5:$mapq $6:$cigar $7:$rnext $8:$pnext $9:$tlen $10:$seq $11:$qual vcf: $1:$chrom $2:$pos $3:$id $4:$ref $5:$alt $6:$qual $7:$filter $8:$info gff: $1:$seqname $2:$source $3:$feature $4:$start $5:$end $6:$score $7:$filter $8:$strand $9:$group $10:$attribute fastx: $1:$name $2:$seq $3:$qual $4:$comment
M120_S2_L001_R1_001_fastqc folder with an html page fastqc_report.html inside.
You can load it up in your browser to assess your data through graphs and summary tables.
[uzi@quince-srv2 ~/linuxTutorial]$ fastqc data/M120_*R1*.fastq
Started analysis of M120_S2_L001_R1_001.fastq
Approx 5% complete for M120_S2_L001_R1_001.fastq
Approx 10% complete for M120_S2_L001_R1_001.fastq
Approx 15% complete for M120_S2_L001_R1_001.fastq
Approx 20% complete for M120_S2_L001_R1_001.fastq
Approx 25% complete for M120_S2_L001_R1_001.fastq
Approx 30% complete for M120_S2_L001_R1_001.fastq
Approx 35% complete for M120_S2_L001_R1_001.fastq
Approx 40% complete for M120_S2_L001_R1_001.fastq
Approx 45% complete for M120_S2_L001_R1_001.fastq
Approx 50% complete for M120_S2_L001_R1_001.fastq
Approx 55% complete for M120_S2_L001_R1_001.fastq
Approx 60% complete for M120_S2_L001_R1_001.fastq
Approx 65% complete for M120_S2_L001_R1_001.fastq
Approx 70% complete for M120_S2_L001_R1_001.fastq
Approx 75% complete for M120_S2_L001_R1_001.fastq
Approx 80% complete for M120_S2_L001_R1_001.fastq
Approx 85% complete for M120_S2_L001_R1_001.fastq
Approx 90% complete for M120_S2_L001_R1_001.fastq
Approx 95% complete for M120_S2_L001_R1_001.fastq
Approx 100% complete for M120_S2_L001_R1_001.fastq
Analysis complete for M120_S2_L001_R1_001.fastq
[uzi@quince-srv2 ~/linuxTutorial]$
For example, here is the file generated for the above M120_S2_L001_R1_001.fastq file:
[uzi@quince-srv2 ~/linuxTutorial/data]$ bwa index Cdiff078.fa
Use BWA-MEM to align paired-end sequences. Briefly, the algorithm works by seeding alignments with maximal exact matches (MEMs) and then extending seeds with the affine-gap Smith-Waterman algorithm (SW).
From BWA doc, it is suggested that for 70bp or longer Illumina, 454, Ion Torrent and Sanger reads, assembly contigs and BAC sequences, BWA-MEM is usually the preferred algorithm. For short sequences, BWA-backtrack may be better. BWA-SW may have better sensitivity when alignment gaps are frequent.
[uzi@quince-srv2 ~/linuxTutorial]$ bwa mem data/Cdiff078.fa data/M120_*R1*.fastq data/M120_*R2*.fastq > aln-pe.sam
We have generated a sam file (aln-pe.sam) which consist of two types of lines: headers and alignments. Headers begin with @, and provide meta-data regarding the entire alignment file.
Alignments begin with any character except @, and describe a single alignment of a sequence read against the reference genome.
Note that each read in a FASTQ file may align to multiple regions within a reference genome, and an individual read can therefore result in multiple alignments.
In the SAM format, each of these alignments is reported on a separate line. Also, each alignment has 11 mandatory fields, followed by a variable number of optional fields. Each of the fields is described in the table below:
| Col | Field | Description |
| 1 | QNAME | Query template/pair NAME |
| 2 | FLAG | bitwise FLAG |
| 3 | RNAME | Reference sequence NAME |
| 4 | POS | 1-based leftmost POSition/coordinate of clipped sequence |
| 5 | MAPQ | MAPping Quality (Phred-scaled) |
| 6 | CIAGR | extended CIGAR string |
| 7 | MRNM | Mate Reference sequence NaMe (‘=’ if same as RNAME) |
| 8 | MPOS | 1-based Mate POSistion |
| 9 | TLEN | inferred Template LENgth (insert size) |
| 10 | SEQ | query SEQuence on the same strand as the reference |
| 11 | QUAL | query QUALity (ASCII-33 gives the Phred base quality) |
| 12+ | OPT | variable OPTional fields in the format TAG:VTYPE:VALUE |
| Flag | Chr | Description |
| 0x0001 | p | the read is paired in sequencing |
| 0x0002 | P | the read is mapped in a proper pair |
| 0x0004 | u | the query sequence itself is unmapped |
| 0x0008 | U | the mate is unmapped |
| 0x0010 | r | strand of the query (1 for reverse) |
| 0x0020 | R | strand of the mate |
| 0x0040 | 1 | the read is the first read in a pair |
| 0x0080 | 2 | the read is the second read in a pair |
| 0x0100 | s | the alignment is not primary |
| 0x0200 | f | the read fails platform/vendor quality checks |
| 0x0400 | d | the read is either a PCR or an optical duplicate |
aln-pe.sam file:
[uzi@quince-srv2 ~/linuxTutorial]$ samstat aln-pe.sam
Running the above code will generate a aln-pe.sam.samstat.html file which you can open in your browser (be patient, it takes a bit of time to load). Plots such as "Reads Length Distributions" and "Base Quality Distributions" may be of interest to you:
[uzi@quince-srv2 ~/linuxTutorial]$ samtools view -h -b -S aln-pe.sam > aln-pe.bam
Extract only those sequences that were mapped against the reference database. Use -F 4 switch.
[uzi@quince-srv2 ~/linuxTutorial]$ samtools view -b -F 4 aln-pe.bam > aln-pe.mapped.bam
Generate a file lengths.genome that contains two entries per row: genome identifier and the corresponding genome length:
[uzi@quince-srv2 ~/linuxTutorial]$ samtools view -H aln-pe.mapped.bam | perl -ne 'if ($_ =~ m/^\@SQ/) { print $_ }' | perl -ne 'if ($_ =~ m/SN:(.+)\s+LN:(\d+)/) { print $1, "\t", $2, "\n"}' > lengths.genome [uzi@quince-srv2 ~/linuxTutorial]$ cat lengths.genome Cdiff078_C01 9165 Cdiff078_C02 93786 Cdiff078_C03 752 Cdiff078_C04 5361 Cdiff078_C05 70058 Cdiff078_C06 23538 Cdiff078_C07 98418 Cdiff078_C08 361074 Cdiff078_C09 45183 Cdiff078_C10 141523 Cdiff078_C11 21992 Cdiff078_C12 2353 Cdiff078_C13 133975 Cdiff078_C14 3374 Cdiff078_C15 9744 Cdiff078_C16 25480 Cdiff078_C17 293596 Cdiff078_C18 7057 Cdiff078_C19 73989 Cdiff078_C20 248092 Cdiff078_C21 41937 Cdiff078_C22 65693 Cdiff078_C23 21321 Cdiff078_C24 440055 Cdiff078_C25 210910 Cdiff078_C26 164162 Cdiff078_C27 22782 Cdiff078_C28 201701 Cdiff078_C29 13447 Cdiff078_C30 101704 Cdiff078_C31 146436 Cdiff078_C32 61153 Cdiff078_C33 59640 Cdiff078_C34 193273 Cdiff078_C35 18395 Cdiff078_C36 25573 Cdiff078_C37 61616 Cdiff078_C38 4117 Cdiff078_C39 110461 Cdiff078_C40 125351 Cdiff078_C41 38508 Cdiff078_C42 113221 Cdiff078_C43 500 Cdiff078_C44 547 Cdiff078_C45 613 Cdiff078_C46 649 Cdiff078_C47 666 Cdiff078_C48 783 Cdiff078_C49 872 Cdiff078_C50 872 Cdiff078_C51 879 Cdiff078_C52 921 Cdiff078_C53 955 Cdiff078_C54 1217 Cdiff078_C55 1337 Cdiff078_C56 1445 Cdiff078_C57 2081 Cdiff078_C58 2098 Cdiff078_C59 2512 Cdiff078_C60 2800 Cdiff078_C61 4372 [uzi@quince-srv2 ~/linuxTutorial]$Sort BAM file. Many of the downstream analysis programs that use BAM files actually require a sorted BAM file.
-m specifies the maximum memory to use, and can be changed to fit your system.
[uzi@quince-srv2 ~/linuxTutorial]$ samtools sort -m 1000000000 aln-pe.mapped.bam aln-pe.mapped.sorted
We will now use bedtools. It is a very useful suite of programs for working with SAM/BAM, BED, VCF and GFF files, files that you will encouter many times doing NGS analysis. -ibam switch takes indexed bam file that we generated earlier, -d reports the depth at each genome position with 1-based coordinates, and -g is used to provide the genome lengths file we generated earlier.
The coverage flags are explained pictorially from genomecov man page:
http://bedtools.readthedocs.org/en/latest/_images/genomecov-glyph.png
[uzi@quince-srv2 ~/linuxTutorial]$ bedtools genomecov -ibam aln-pe.mapped.sorted.bam -d -g lengths.genome > aln-pe.mapped.bam.perbase.cov
Look at the first few entries in the file generated above. First column is genome identifier, second column is position on genome, and third column is coverage.
[uzi@quince-srv2 ~/linuxTutorial]$ head aln-pe.mapped.bam.perbase.cov Cdiff078_C01 1 41 Cdiff078_C01 2 41 Cdiff078_C01 3 42 Cdiff078_C01 4 42 Cdiff078_C01 5 42 Cdiff078_C01 6 44 Cdiff078_C01 7 44 Cdiff078_C01 8 44 Cdiff078_C01 9 44 Cdiff078_C01 10 44 [uzi@quince-srv2 ~/linuxTutorial]$Now we will count only those positions where we have >0 coverage.
[uzi@quince-srv2 ~/linuxTutorial]$ awk -F"\t" '$3>0{print $1}' aln-pe.mapped.bam.perbase.cov | sort | uniq -c > aln-pe.mapped.bam.perbase.count
To see what we have done, use the cat command
[uzi@quince-srv2 ~/linuxTutorial]$ cat aln-pe.mapped.bam.perbase.count 9165 Cdiff078_C01 93786 Cdiff078_C02 752 Cdiff078_C03 5361 Cdiff078_C04 70058 Cdiff078_C05 23538 Cdiff078_C06 98418 Cdiff078_C07 333224 Cdiff078_C08 44803 Cdiff078_C09 141523 Cdiff078_C10 21969 Cdiff078_C11 2292 Cdiff078_C12 133974 Cdiff078_C13 1762 Cdiff078_C14 50 Cdiff078_C15 10232 Cdiff078_C16 293440 Cdiff078_C17 7057 Cdiff078_C18 73989 Cdiff078_C19 248092 Cdiff078_C20 41937 Cdiff078_C21 65447 Cdiff078_C22 21321 Cdiff078_C23 439123 Cdiff078_C24 210910 Cdiff078_C25 164162 Cdiff078_C26 22782 Cdiff078_C27 201701 Cdiff078_C28 13447 Cdiff078_C29 98510 Cdiff078_C30 146261 Cdiff078_C31 61153 Cdiff078_C32 44523 Cdiff078_C33 193180 Cdiff078_C34 18395 Cdiff078_C35 25573 Cdiff078_C36 61616 Cdiff078_C37 4117 Cdiff078_C38 62897 Cdiff078_C39 125351 Cdiff078_C40 38508 Cdiff078_C41 113221 Cdiff078_C42 442 Cdiff078_C43 649 Cdiff078_C46 663 Cdiff078_C47 766 Cdiff078_C48 580 Cdiff078_C51 1110 Cdiff078_C54 1445 Cdiff078_C56 2512 Cdiff078_C59 2800 Cdiff078_C60 [uzi@quince-srv2 ~/linuxTutorial]$We will now use the above file with
lengths.genome to calculate the proportions of genomes/contigs covered using the following one-liner. It reads lengths.genome line by line, assigns the genome identifier to myArray[0], it's length to myArray[1]. It then searches the identifier in aln-pe.mapped.bam.perbase.count, extracts the base count, and uses bc to calculate the proportions.
[uzi@quince-srv2 ~/linuxTutorial]$ while IFS=$'\t' read -r -a myArray; do echo -e "${myArray[0]},$( echo "scale=5;0"$(awk -v pattern="${myArray[0]}" '$2==pattern{print $1}' aln-pe.mapped.bam.perbase.count)"/"${myArray[1]} | bc ) "; done < lengths.genome > aln-pe.mapped.bam.genomeproportion [uzi@quince-srv2 ~/linuxTutorial]$ cat aln-pe.mapped.bam.genomeproportion Cdiff078_C01,1.00000 Cdiff078_C02,1.00000 Cdiff078_C03,1.00000 Cdiff078_C04,1.00000 Cdiff078_C05,1.00000 Cdiff078_C06,1.00000 Cdiff078_C07,1.00000 Cdiff078_C08,.92286 Cdiff078_C09,.99158 Cdiff078_C10,1.00000 Cdiff078_C11,.99895 Cdiff078_C12,.97407 Cdiff078_C13,.99999 Cdiff078_C14,.52222 Cdiff078_C15,.00513 Cdiff078_C16,.40156 Cdiff078_C17,.99946 Cdiff078_C18,1.00000 Cdiff078_C19,1.00000 Cdiff078_C20,1.00000 Cdiff078_C21,1.00000 Cdiff078_C22,.99625 Cdiff078_C23,1.00000 Cdiff078_C24,.99788 Cdiff078_C25,1.00000 Cdiff078_C26,1.00000 Cdiff078_C27,1.00000 Cdiff078_C28,1.00000 Cdiff078_C29,1.00000 Cdiff078_C30,.96859 Cdiff078_C31,.99880 Cdiff078_C32,1.00000 Cdiff078_C33,.74652 Cdiff078_C34,.99951 Cdiff078_C35,1.00000 Cdiff078_C36,1.00000 Cdiff078_C37,1.00000 Cdiff078_C38,1.00000 Cdiff078_C39,.56940 Cdiff078_C40,1.00000 Cdiff078_C41,1.00000 Cdiff078_C42,1.00000 Cdiff078_C43,.88400 Cdiff078_C44,0 Cdiff078_C45,0 Cdiff078_C46,1.00000 Cdiff078_C47,.99549 Cdiff078_C48,.97828 Cdiff078_C49,0 Cdiff078_C50,0 Cdiff078_C51,.65984 Cdiff078_C52,0 Cdiff078_C53,0 Cdiff078_C54,.91207 Cdiff078_C55,0 Cdiff078_C56,1.00000 Cdiff078_C57,0 Cdiff078_C58,0 Cdiff078_C59,1.00000 Cdiff078_C60,1.00000 Cdiff078_C61,0We have a total of 61 genomes/contigs in the reference database. To see how many genomes/contigs we recovered, we will use the following one-liner:
[uzi@quince-srv2 ~/linuxTutorial]$ awk -F "," '{sum+=$NF} END{print "Total genomes covered:"sum}' aln-pe.mapped.bam.genomeproportion Total genomes covered:47.5224We also need genome/contig coverage, which we can calculate as:
[uzi@quince-srv2 ~/linuxTutorial]$ bedtools genomecov -ibam aln-pe.mapped.sorted.bam -g lengths.genome | awk -F"\t" '!/^genome/{l[$1]=l[$1]+($2 *$3);r[$1]=$4} END {for (i in l){print i","(l[i]/r[i])}}' > aln-pe.mapped.bam.genomecoverage [uzi@quince-srv2 ~/linuxTutorial]$ cat aln-pe.mapped.bam.genomecoverage Cdiff078_C10,61.5467 Cdiff078_C11,68.9158 Cdiff078_C12,79.7875 Cdiff078_C13,61.2645 Cdiff078_C14,57.3438 Cdiff078_C15,0.0812808 Cdiff078_C16,23.5227 Cdiff078_C17,57.358 Cdiff078_C30,59.3333 Cdiff078_C18,55.5597 Cdiff078_C31,62.147 Cdiff078_C19,56.3139 Cdiff078_C32,66.0493 Cdiff078_C33,48.8165 Cdiff078_C34,65.7106 Cdiff078_C35,62.7728 Cdiff078_C36,62.7535 Cdiff078_C37,67.2169 Cdiff078_C51,1.05916 Cdiff078_C38,61.9871 Cdiff078_C39,37.3289 Cdiff078_C54,6.46754 Cdiff078_C56,815.224 Cdiff078_C59,801.998 Cdiff078_C01,67.3333 Cdiff078_C02,67.4621 Cdiff078_C03,103.848 Cdiff078_C04,65.4128 Cdiff078_C05,66.1244 Cdiff078_C06,66.239 Cdiff078_C07,76.0081 Cdiff078_C20,55.6661 Cdiff078_C08,60.4236 Cdiff078_C21,56.2321 Cdiff078_C09,76.9986 Cdiff078_C22,56.8815 Cdiff078_C23,53.2772 Cdiff078_C24,56.9991 Cdiff078_C25,57.4446 Cdiff078_C26,59.296 Cdiff078_C40,66.0074 Cdiff078_C27,59.4391 Cdiff078_C41,67.5941 Cdiff078_C28,59.8319 Cdiff078_C42,69.4415 Cdiff078_C29,60.961 Cdiff078_C43,4.812 Cdiff078_C46,29.3837 Cdiff078_C60,62.1336 Cdiff078_C47,7.95946 Cdiff078_C48,15.3436 [uzi@quince-srv2 ~/linuxTutorial]$Sort the original bam file
[uzi@quince-srv2 ~/linuxTutorial]$ samtools sort -m 1000000000 aln-pe.bam aln-pe.sorted
Now we will check alignment statistics using the Picard tools. Note that the awk statement given below is used to transpose the original table and you can do without it.
[uzi@quince-srv2 ~/linuxTutorial]$ java -jar $(which CollectAlignmentSummaryMetrics.jar) INPUT=aln-pe.sorted.bam OUTPUT=aln-pe.sorted.alignment_stats.txt REFERENCE_SEQUENCE=data/Cdiff078.fa [uzi@quince-srv2 ~/linuxTutorial]$ grep -vi -e "^#" -e "^$" aln-pe.sorted.alignment_stats.txt | awk -F"\t" '{ for (i=1; i<=NF; i++) {a[NR,i] = $i}}NF>p{p=NF}END{for(j=1;j<=p;j++){str=a[1,j];for(i=2; i<=NR; i++){str=str"\t"a[i,j];} print str}}' CATEGORY FIRST_OF_PAIR SECOND_OF_PAIR PAIR TOTAL_READS 425271 425038 850309 PF_READS 425271 425038 850309 PCT_PF_READS 1 1 1 PF_NOISE_READS 0 0 0 PF_READS_ALIGNED 407011 405258 812269 PCT_PF_READS_ALIGNED 0.957063 0.953463 0.955263 PF_ALIGNED_BASES 119451610 118113100 237564710 PF_HQ_ALIGNED_READS 401018 399295 800313 PF_HQ_ALIGNED_BASES 118606615 117274833 235881448 PF_HQ_ALIGNED_Q20_BASES 116971078 111640501 228611579 PF_HQ_MEDIAN_MISMATCHES 0 0 0 PF_MISMATCH_RATE 0.002359 0.007186 0.004759 PF_HQ_ERROR_RATE 0.002269 0.007065 0.004653 PF_INDEL_RATE 0.000124 0.00013 0.000127 MEAN_READ_LENGTH 299.093366 298.832657 298.963048 READS_ALIGNED_IN_PAIRS 404714 404545 809259 PCT_READS_ALIGNED_IN_PAIRS 0.994356 0.998241 0.996294 BAD_CYCLES 0 0 0 STRAND_BALANCE 0.500072 0.500484 0.500278 PCT_CHIMERAS 0.014823 0.014668 0.014746 PCT_ADAPTER 0.000285 0.000261 0.000273 SAMPLE LIBRARY READ_GROUP [uzi@quince-srv2 ~/linuxTutorial]$The detailed description of these summary metrics are given here. From this link,
PF_MISMATCH_RATE, PF_HQ_ERROR_RATE, and PF_INDEL_RATE are of interest to us.
As can be seen, the error rates are quite low and we can proceed with the analysis. Next we would like to calculate GC bias. For this purpose, we will index aln-pe.mapped.sorted.bam file.
[uzi@quince-srv2 ~/linuxTutorial]$ samtools index aln-pe.mapped.sorted.bam [uzi@quince-srv2 ~/linuxTutorial]$ for i in $(samtools view -H aln-pe.mapped.sorted.bam | awk -F"\t" '/@SQ/{gsub("^SN:","",$2);print $2}' );do samtools view -b aln-pe.mapped.sorted.bam $i > aln-pe.mapped.sorted.$i.bam; java -Xmx2g -jar $(which CollectGcBiasMetrics.jar) R=data/Cdiff078.fa I=aln-pe.mapped.sorted.$i.bam O=aln-pe.mapped.sorted.${i}_GCBias.txt CHART=aln-pe.mapped.sorted.${i}_GCBias.pdf ASSUME_SORTED=true; doneIn the above one-liner,
CollectGcBiasMetrics.jar will generate a GC bias plot for each contig, and will look like these:
[uzi@quince-srv2 ~/linuxTutorial]$ for i in $(ls *_GCBias.txt); do awk -v k="$i" '!/^#/ && !/^$/ && !/^GC/ && !/?/{print k"\t"$1"\t"$5}' $i; done | perl -ane '$r{$F[0].":".$F[1]}=$F[2];unless($F[0]~~@s){push @s,$F[0];}unless($F[1]~~@m){push @m,$F[1];}END{print "Contigs\t".join("\t",@s)."\n";for($i=0;$i<@m;$i++){print $m[$i];for($j=0;$j<@s;$j++){(not defined $r{$s[$j].":".$m[$i]})?print "\t".0:print"\t".$r{$s[$j].":".$m[$i]};}print "\n";}}' | sed '1s/aln-pe\.mapped\.sorted\.//g;1s/_GCBias\.txt//g' > aln-pe.mapped.sorted.bam.gcbias
GENERATEtable.sh script. If the data stream is of the form [Contig]\t[Feature]\t[Value], then you can pipe the stream to GENERATEtable.sh to obtain a Contig X Feature table:
$ cat test.tsv contig1 F1 12.2 contig1 F2 34.2 contig1 F3 45.2 contig2 F2 56.3 contig2 F3 56.2 contig3 F1 45.4 contig3 F2 56.3 contig4 F1 23.5 contig5 F1 24.5 $ cat GENERATEtable.sh #!/bin/bash less <&0| \ perl -ane '$r{$F[0].":".$F[1]}=$F[2]; unless($F[0]~~@s){ push @s,$F[0];} unless($F[1]~~@m){ push @m,$F[1];} END{ print "Contigs\t".join("\t",@s)."\n"; for($i=0;$i<@m;$i++){ print $m[$i]; for($j=0;$j<@s;$j++){ (not defined $r{$s[$j].":".$m[$i]})?print "\t".0:print"\t".$r{$s[$j].":".$m[$i]};} print "\n";}}' $ cat test.tsv | ./GENERATEtable.sh Contigs contig1 contig2 contig3 contig4 contig5 F1 12.2 0 45.4 23.5 24.5 F2 34.2 56.3 56.3 0 0 F3 45.2 56.2 0 0 0
[uzi@quince-srv2 ~/linuxTutorial]$ head aln-pe.mapped.sorted.bam.gcbias Contigs Cdiff078_C01 Cdiff078_C02 Cdiff078_C03 Cdiff078_C04 Cdiff078_C05 Cdiff078_C06 Cdiff078_C07 Cdiff078_C08 Cdiff078_C09 Cdiff078_C10 Cdiff078_C11 Cdiff078_C12 Cdiff078_C13 Cdiff078_C14 Cdiff078_C15 Cdiff078_C16 Cdiff078_C17 Cdiff078_C18 Cdiff078_C19 Cdiff078_C20 Cdiff078_C21 Cdiff078_C22 Cdiff078_C23 Cdiff078_C24 Cdiff078_C25 Cdiff078_C26 Cdiff078_C27 Cdiff078_C28 Cdiff078_C29 Cdiff078_C30 Cdiff078_C31 Cdiff078_C32 Cdiff078_C33 Cdiff078_C34 Cdiff078_C35 Cdiff078_C36 Cdiff078_C37 Cdiff078_C38 Cdiff078_C39 Cdiff078_C40 Cdiff078_C41 Cdiff078_C42 Cdiff078_C43 Cdiff078_C46 Cdiff078_C47 Cdiff078_C48 Cdiff078_C51 Cdiff078_C54 Cdiff078_C56 Cdiff078_C59 Cdiff078_C60 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 19.855622 0 113.514035 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 5.94012 0 0 9.265957 0 21.189287 0 0 1.013634 0 0 0 0 0 0 0 0 0 0 21.79477 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0 0 2.099257 0 0 8.004623 17.472841 0 0 4.059842 0.50366 0 0 1.011679 0 0 3.954422 0 15.404681 0 4.36693 0 1.355168 0 0 0 0 0 0 0 0 0 0 7 0 0 0 0 1.898603 0 0 0 5.621592 3.322151 0 0 1.063016 0 0 0 0.395449 0 0 1.454454 11.850605 1.167889 0 2.616454 0.919685 0.229238 0 0.554199 0 0 2.647623 1.10446 0 0.699548 5.625801 0 3.721214 0 1.484725 0 0 0 0 0 0 0 0 0 23.930557 8 0 0.647487 0 0 5.277073 0 1.787893 0 4.918962 1.302196 5.292116 0 0.73865 2.481033 0 0 0.488503 0 0 1.796704 8.023401 0.270507 0 2.181692 1.27811 1.168119 0 1.625944 0 0 1.672486 1.534896 0.32406 4.34351 0 0.738781 0 0.917049 0 0 0.646285 0 0 0 0 0 9 0 0.2347 0 0 1.275217 0 0.465283 0.362786 3.35627 0.457715 2.237992 0 1.725467 4.796381 0 0 0.885359 0 0.719303 1.9538 3.673656 0.915163 1.311588 1.97216 0.494173 1.283086 0 0.909904 0 0.080909 1.131649 2.225468 0 0.352394 1.259544 0 1.071169 13.719311 1.329643 0.540768 0.765043 0.562234 0 0 0 0 0 0 0 0 8.036611 10 0.818789 0.522956 0 0 1.092856 0.491785 0.512673 0.769493 1.150524 0.794326 2.959123 0 1.437926 1.849716 0 0 0.713225 0.63983 0.739728 1.748817 2.30876 1.411724 0.674416 1.696827 1.355215 1.504252 0.187891 1.361078 0 0.582441 0.997528 2.288661 0 1.147602 0 1.42304 1.223984 4.409023 0.569749 0.463436 0.393384 0.321222 0 0 0 0 0 0 0 0 11 2.714431 0.315217 0 0 0.724603 0.494047 0.247215 0.885474 1.083578 0.904377 2.378187 0 1.954739 0 0 0.50476 0.713456 0 1.226165 1.513891 1.138604 2.795913 1.355035 1.755064 1.697555 1.208937 0.566266 1.204967 0 0.927834 1.018818 1.456152 0 1.04366 0.260253 0.857751 0.848433 0.892897 0.335209 0.829905 0.755116 0 0 0 0 0 0 0 0 2.490848You can then use the following
R code to generate a GC vs Coverage table which shows that at very GC, coverages go down (note that these are the smoothed values across all genomes/contigs):
library(ggplot2) library(reshape) data_table <- read.csv("aln-pe.mapped.sorted.bam.gcbias",header=TRUE,row.names=1,sep="\t") df<-NULL for(i in names(data_table)){ tmp<-data.frame(rownames(data_table),data_table[,i],rep(i,dim(data_table)[1])) if(is.null(df)){df<-tmp}else{df<-rbind(df,tmp)} } names(df)<-c("GC","Coverage","Contigs") df$GC<-as.numeric(df$GC) p<-ggplot(df,aes(GC,Coverage,group=Contigs)) + geom_smooth(aes(group="dummy"),method = "loess", formula = y ~ x, size = 1)+ theme_bw()+ theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.5))+ pdf("aln-pe.mapped.sorted.bam.gcbias.pdf") print(p) dev.off()
[uzi@quince-srv2 ~/linuxTutorial]$ java -Xmx2g -jar $(which MeanQualityByCycle.jar) INPUT=aln-pe.mapped.sorted.bam OUTPUT=aln-pe.mapped.sorted.mqc.txt CHART_OUTPUT=aln-pe.mapped.sorted.mqc.pdf
[uzi@quince-srv2 ~/linuxTutorial]$ java -Xmx2g -jar $(which QualityScoreDistribution.jar) INPUT=aln-pe.mapped.sorted.bam OUTPUT=aln-pe.mapped.sorted.qsd.txt CHART_OUTPUT=aln-pe.mapped.sorted.qsd.pdf
input.txt file which contains listing of BAM files we want to compare. To save time, I am only considering 3 files:
[uzi@quince-srv2 ~/linuxTutorial]$ awk '{split($0,k,".");print k[4]"\t"$0}' <(ls *.sorted.Cdiff078_C4*.bam | head -3) > input.txt [uzi@quince-srv2 ~/linuxTutorial]$ cat input.txt Cdiff078_C40 aln-pe.mapped.sorted.Cdiff078_C40.bam Cdiff078_C41 aln-pe.mapped.sorted.Cdiff078_C41.bam Cdiff078_C42 aln-pe.mapped.sorted.Cdiff078_C42.bam [uzi@quince-srv2 ~/linuxTutorial]$ qualimap multi-bamqc -d input.txt -r -outdir qualimapThis will generate a
qualimap folder that will contain a multisampleBamQcReport.html file that you can load in your browser. Do check the "PCA" and "Insert Size Histogram" plot.
samtools flagstat to get mapping statistics:
[uzi@quince-srv2 ~/linuxTutorial]$ samtools flagstat aln-pe.bam 850309 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 duplicates 812269 + 0 mapped (95.53%:-nan%) 850309 + 0 paired in sequencing 425271 + 0 read1 425038 + 0 read2 795935 + 0 properly paired (93.61%:-nan%) 809259 + 0 with itself and mate mapped 3010 + 0 singletons (0.35%:-nan%) 11922 + 0 with mate mapped to a different chr 8256 + 0 with mate mapped to a different chr (mapQ>=5)With
Bioawk, you can do amazing things: