Six weeks basic computing workshop


Week 4: samtools, bowtie2, bedtools multicov, igv

($ signs at the begining of the command lines indicate after $ is a command)
  1. Connect to the cluster using ssh.
  2. $ ssh your_user@ghpcc06.umassrc.org 
  3. Prepare your working directory.
  4. $ clear 
    $ mkdir ~/class/week4
    $ cd ~/class/week4
  5. Create necessary directories
  6. $ mkdir seq 
    $ mkdir index 
    $ mkdir annot 
  7. Copy example files to their directories
  8. $ cp /share/data/umw_biocore/genome_data/mousetest/mm10/example/seq/* seq/. 
    $ cp /share/data/umw_biocore/genome_data/mousetest/mm10/example/annot/* annot/. 
    $ cp /share/data/umw_biocore/genome_data/mousetest/mm10/example/index/* index/. 
  9. First run fastqc for quality checks of the libraries
  10. $ mkdir fastqc_results
    $ module load fastqc/0.10.1
    $ fastqc seq/example_qc.fastq -o fastqc_results
  11. Now lets build an index file for mm10 example genome
  12. $ module load bowtie2/2-2.1.0
    $ bowtie2-build index/mm10.fa index/mm10
  13. We can now align our sequences to the reference genome
  14. mkdir mappings
    $ bowtie2 -N 1 -k 1 --no-unal -x index/mm10 -1 seq/control_rep1_R1.fastq -2 seq/control_rep1_R2.fastq > mappings/control_rep1.sam
  15. To be able to find read counts per gene, we need to convert sam file to bam file and need sort and index. Since bedtools requires sorted and indexed files to be able to find gene counts.
  16.  $ module load samtools/0.0.19 
     $ samtools view -bT index/mm10.fa mappings/control_rep1.sam > mappings/control_rep1.bam 
     $ samtools sort mappings/control_rep1.bam mappings/control_rep1.sorted
     $ samtools index mappings/control_rep1.sorted.bam
  17. Now we can count number of reads per genes
  18. $ bedtools multicov -p -F -bams mappings/control_rep1.sorted.bam -bed annot/mm10ref.bed | awk '{print $4"\t"$13}' > readcounts.txt 
  19. To visualize a bam file in igv please downlaod igv program
  20. Mac Users: IGV_2.3.34.app.zip

    Windows/Linux Users: IGV_2.3.34.zip

  21. Download bam file and index file to your computer
  22. Start igv and select mouse(mm10)
  23. Load the bam file
  24. Your homework (deadline: Aug 6, Thursday 17:00pm):


    1. Please create “homework/week4” directory in your home folder.

    2. Download all paired end sequence files from /share/data/umw_biocore/genome_data/mousetest/mm10/example/seq/.

    3. Map them those sequence files to mm10 genome sequence file is located in the following directory /share/data/umw_biocore/genome_data/mousetest/mm10/example/index

    4. Count the number of reads per isoform and report them in a 8 column count file that includes following columns; gene_id gene_len control_rep1 control_rep2 ...
    5. Write all the command/s you used into Readme.txt file in the same folder.