RNA-Seq File Format Practicals

A primer to the different types of RNA-Seq data files

Before you start

1. Check the github page for this document.
2. If you would like to follow this practical from your own machine, you may need to install some of the software tools (follow relevant links) or linux commands.
3. If 2 is true, you may also need to set up your $PATH environment variable

Aims

1. To prepare other classes ahead
2. To gain confidence

• of processing RNA-Seq data files
• or linux (or unix) command lines

You will learn

1. Various file formats

• FASTA and FASTQ
• GTF and or GFF
• SAM and BAM
• BED

2. Accessing public or published data

• Ensembl
• NCBI SRA

3. And related linux commands and software tools

• awk
• sort
• zcat or zless
• samtools
• bedtools
• sra toolkits

You will NOT learn (from me)

1. Various RNA-Seq data processing and analysis, such as

• Quality controls of sequenced reads
• Adaptor (and/or poor quality read) trimming
• Alignment (mapping) of reads
• Quantification of transcript-level (or gene-level) abundance
• Differentially expressed gene/transcript/exon analysis
• Transcriptome reconstruction (or assembly)

2. Detailed use of software tools related with above

• fastqc
• cutadapt
• TopHat2 and HiSat2
• featureCounts, Salmon Kallisto
• DESeq2, DEXSeq, and edgeR
• StringTie and TACO

File Formats

https://genome.ucsc.edu/FAQ/FAQformat

FASTA

• https://en.wikipedia.org/wiki/FASTA_format

• The letters should follow the IUPAC codes

• Task 1: download the FASTA file of human chr22 (from here) and calculate the size (i.e. length) of the chromosome

wget ftp://ftp.ensembl.org/pub/release-95/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.22.fa.gz

zless Homo_sapiens.GRCh38.dna.chromosome.22.fa.gz | head -n5

>22 dna:chromosome chromosome:GRCh38:22:1:50818468:1 REF
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

zcat Homo_sapiens.GRCh38.dna.chromosome.22.fa.gz | awk '!/^>/{LEN+=length($0)}END{print LEN}' 

50818468

• Task 2: get the sequence from 40000001 to 40000300 of chr22

zcat Homo_sapiens.GRCh38.dna.chromosome.22.fa.gz | awk 'BEGIN{min=50818468;max=0;start=40000001;end=40000300}!/^>/{LEN+=length($0); if(LEN>=start && LEN-length($0)+1<=end){min=(min>LEN-length($0)+1)?LEN-length($0)+1:min; max=(max>LEN)?max:LEN; SEQ=SEQ$0; print LEN-length($0)+1,LEN,$0,length($0)} }END{print substr(SEQ,start-min+1,end-start+1)}'

39999961 40000020 TTGGGAGAGTGCATTTGACCCTCCTGGATTGGGAACTTTGGAAAAGAGAGACTGGCTTAC 60
40000021 40000080 TCTCTCAGTTTCCTTTCCGGTACCCGAACAGTGCTCTGACTGTCACAGGTGATCAATTAT 60
40000081 40000140 CTGTTTATAGAATTGAAAACTGTCTTCTATTGGCCGCCTGTCCCTCTGGGGAAGCAAGTG 60
40000141 40000200 TCTGGAGATGTCCCCATGGAATTGCACATTCCCTGGGGAACTCCTTTAATCTACATCGTA 60
40000201 40000260 ACAAGCTGGCTGGCAGCAGAGAGAAGAATGCTCTAAATTGGGGATATGGGGTGCAGGGAG 60
40000261 40000320 GGAGGTGCATCAGGGGATTTAGGGAAAGGCTCCAGAGCAGGAGGCATGGCTGGCTGAGGT 60
GAAAAGAGAGACTGGCTTACTCTCTCAGTTTCCTTTCCGGTACCCGAACAGTGCTCTGACTGTCACAGGTGATCAATTATCTGTTTATAGAATTGAAAACTGTCTTCTATTGGCCGCCTGTCCCTCTGGGGAAGCAAGTGTCTGGAGATGTCCCCATGGAATTGCACATTCCCTGGGGAACTCCTTTAATCTACATCGTAACAAGCTGGCTGGCAGCAGAGAGAAGAATGCTCTAAATTGGGGATATGGGGTGCAGGGAGGGAGGTGCATCAGGGGATTTAGGGAAAGGCTCCAGAGCAG

• Compare the result above and the following method (using the ensembl rest api)

curl https://rest.ensembl.org/sequence/region/human/22:40000001..40000300:1?content-type=text/plain

  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed

  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0
100   300  100   300    0     0   2882      0 --:--:-- --:--:-- --:--:--  2912
GAAAAGAGAGACTGGCTTACTCTCTCAGTTTCCTTTCCGGTACCCGAACAGTGCTCTGACTGTCACAGGTGATCAATTATCTGTTTATAGAATTGAAAACTGTCTTCTATTGGCCGCCTGTCCCTCTGGGGAAGCAAGTGTCTGGAGATGTCCCCATGGAATTGCACATTCCCTGGGGAACTCCTTTAATCTACATCGTAACAAGCTGGCTGGCAGCAGAGAGAAGAATGCTCTAAATTGGGGATATGGGGTGCAGGGAGGGAGGTGCATCAGGGGATTTAGGGAAAGGCTCCAGAGCAG

• Task 3: what is the sequence of 40000001 to 40000100 of chr21 (using both awk and curl)?

curl https://rest.ensembl.org/sequence/region/human/21:40000001..40000100:1?content-type=text/plain

  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed

  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0
100   100  100   100    0     0   1209      0 --:--:-- --:--:-- --:--:--  1219
GTCTTCCTGACTTCCAAGGTTGCCATGAATCCTTGGTGTTCCTTGGCTTGCAGATGCATTGCCATCTTCCCTGCTTGTCAGTCTGTGTTTCCTTCTTCTT

FASTQ

• https://en.wikipedia.org/wiki/FASTQ_format

• Examples from Gong et al. JCI-I 2018

• Download one of down-sampled (to 1K) FASTQ available from below

wget https://www.dropbox.com/s/h2lzzn94o9zpv26/SLX-9168.D701_D501_1K.C6H3UANXX.s_1.r_1.fq.gz

• Question: how many entries (it’s abvious)

zcat  SLX-9168.D701_D501_1K.C6H3UANXX.s_1.r_1.fq.gz | echo $((`wc -l`/4))

1000

• Let’s look at the first entry (i.e. the first 4 lines) and check the meta information

zcat SLX-9168.D701_D501_1K.C6H3UANXX.s_1.r_1.fq.gz | head -n4

@HISEQ:249:C6H3UANXX:1:1101:1452:2086 1:N:0:NTTACTCGNATAGCCT
NTTTCCAAAAAAATGAAACACACATTAGATATAAACCATCCAAAACTTTAACGATTGAAAACTGTATTCAAAATTAAATGACCCAGAGTAGGTCTCTTTCTGAAAAGTTCCCTTCATATCGGCAA
+
#==AAFGFFEGFDGFCF1<:CF>FGGCEGGGG>G>FBFGGEG@F@FGG1<1FEFGGGG1EGGGGGGGGGGGGCGGFDEGGEE>FD:0FFF::FFGGGGGGGEGGG0ECGGGGDGGG>FGF<8CAG

• Question: sequencing length?

zless SLX-9168.D701_D501_1K.C6H3UANXX.s_1.r_1.fq.gz | head -n 4 | awk 'NR%4==2{print length($0)}'

125

• Check the base quality of the read above from the samformat.info

• Check the base quality using fastqc

fastqc SLX-9168.D701_D501_1K.C6H3UANXX.s_1.r_1.fq.gz

• Another example Ngo et al. Science 2018 (SPR130149)

• Downlaod the first 1K spots from the smallest RNA-Seq fastq files from SRR6481086 (you may need to install the sra toolkit)

fastq-dump -X 1000 --split-files --read-filter pass --outdir ./ --gzip SRR6481086 

zcat SRR6481086_pass_1.fastq.gz | echo $((`wc -l`/4))
zcat SRR6481086_pass_2.fastq.gz | echo $((`wc -l`/4))

1000
1000

• Let’s look at the fist 8 lines

zcat SRR6481086_pass_1.fastq.gz | head -n8

@SRR6481086.1 1 length=75
GGCCATGCAAGATTCCCATTCTTGCGACCCGGGTTCGTTTCCCGGGCGGCGCACCAGATCGGAAGAGCACACGNC
+SRR6481086.1 1 length=75
AAAAAEEEA/EAEEEEE6EEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEE/EEEE6EEEE6AE6EEEEEEE#E
@SRR6481086.2 2 length=75
GAAGAGCACACGTCTGAACTCCAGTCACTAATGCGCATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAAAGGGNG
+SRR6481086.2 2 length=75
AAAA/EE/E6EEEEEE66EEEE/EEE/EEE6EEEAE/A6EEEE6EEEEEE6EEAEE/EEEEEE6EEEEE//AE#E

zcat SRR6481086_pass_2.fastq.gz | head -n8

@SRR6481086.1 1 length=75
GGTGCGCCGCCCNGGCCCCGCNNNCGNGTCGCCCGCCTNNNCNTCTTGCCTGNNCCGCTCGGCCGCGCGTCGTGT
+SRR6481086.1 1 length=75
AAAAAEEEEEEE#EEA//EEE###EE#EEEE/EEEE/E###/#EAEEE//EE##AAEAE/EEA/EEE/EE/E/EA
@SRR6481086.2 2 length=75
GGGGGGGGGGGGNGGGGGGGGNNNGGNGGGGGGGGGGGNNNGNGGGGGGGGGNNGGGGGGGGGGGGGGGGGGGGG
+SRR6481086.2 2 length=75
AAAAAEEEEEEE#EEEEEEEE###EE#EEEEEEEEEEE###E#EEAEAEEEE##EEEEEEEEEEEEEAAAEEAAE

• Check the base quality from the samformat.info

• Check the base quality using fastqc

fastqc SRR6481086_pass_1.fastq.gz SRR6481086_pass_2.fastq.gz

GTF (GFF)

• https://genome.ucsc.edu/FAQ/FAQformat#format3

• https://www.ensembl.org/info/website/upload/gff.html

• Let’s download the Ensembl transcript annotation file

wget ftp://ftp.ensembl.org/pub/release-95/gtf/homo_sapiens/Homo_sapiens.GRCh38.95.gtf.gz

• How many genes and transcripts?

zcat Homo_sapiens.GRCh38.95.gtf.gz | awk '$3=="gene"{GENE++}$3=="transcript"{TRANS++}$3=="exon"{EXON++}END{print GENE,TRANS,EXON}'

58735 206601 1262162

• What else in the file?

zcat Homo_sapiens.GRCh38.95.gtf.gz | awk '!/^#/{print $3}' | sort | uniq -c | sort -k1,1n

    120 Selenocysteine
  58735 gene
  78562 stop_codon
  86454 start_codon
 148491 three_prime_utr
 149930 five_prime_utr
 206601 transcript
 746504 CDS
1262162 exon

• Read the first 10 lines (excluding comments (i.e. ‘#’)) using Bioconductor rtracklayer package

# start your R session using 'R' command from the terminal
rtracklayer::import(format="gtf", text=system("zcat Homo_sapiens.GRCh38.95.gtf.gz | grep -v ^# | head -n 10", intern=T))

Warning in readGFF(filepath, version = version, filter = filter): connection is not positioned at the start of the file, rewinding
  it

GRanges object with 10 ranges and 19 metadata columns:
       seqnames         ranges strand |   source       type     score
          <Rle>      <IRanges>  <Rle> | <factor>   <factor> <numeric>
   [1]        1 [11869, 14409]      + |   havana       gene      <NA>
   [2]        1 [11869, 14409]      + |   havana transcript      <NA>
   [3]        1 [11869, 12227]      + |   havana       exon      <NA>
   [4]        1 [12613, 12721]      + |   havana       exon      <NA>
   [5]        1 [13221, 14409]      + |   havana       exon      <NA>
   [6]        1 [12010, 13670]      + |   havana transcript      <NA>
   [7]        1 [12010, 12057]      + |   havana       exon      <NA>
   [8]        1 [12179, 12227]      + |   havana       exon      <NA>
   [9]        1 [12613, 12697]      + |   havana       exon      <NA>
  [10]        1 [12975, 13052]      + |   havana       exon      <NA>
           phase         gene_id gene_version   gene_name gene_source
       <integer>     <character>  <character> <character> <character>
   [1]      <NA> ENSG00000223972            5     DDX11L1      havana
   [2]      <NA> ENSG00000223972            5     DDX11L1      havana
   [3]      <NA> ENSG00000223972            5     DDX11L1      havana
   [4]      <NA> ENSG00000223972            5     DDX11L1      havana
   [5]      <NA> ENSG00000223972            5     DDX11L1      havana
   [6]      <NA> ENSG00000223972            5     DDX11L1      havana
   [7]      <NA> ENSG00000223972            5     DDX11L1      havana
   [8]      <NA> ENSG00000223972            5     DDX11L1      havana
   [9]      <NA> ENSG00000223972            5     DDX11L1      havana
  [10]      <NA> ENSG00000223972            5     DDX11L1      havana
                             gene_biotype   transcript_id
                              <character>     <character>
   [1] transcribed_unprocessed_pseudogene            <NA>
   [2] transcribed_unprocessed_pseudogene ENST00000456328
   [3] transcribed_unprocessed_pseudogene ENST00000456328
   [4] transcribed_unprocessed_pseudogene ENST00000456328
   [5] transcribed_unprocessed_pseudogene ENST00000456328
   [6] transcribed_unprocessed_pseudogene ENST00000450305
   [7] transcribed_unprocessed_pseudogene ENST00000450305
   [8] transcribed_unprocessed_pseudogene ENST00000450305
   [9] transcribed_unprocessed_pseudogene ENST00000450305
  [10] transcribed_unprocessed_pseudogene ENST00000450305
       transcript_version transcript_name transcript_source
              <character>     <character>       <character>
   [1]               <NA>            <NA>              <NA>
   [2]                  2     DDX11L1-202            havana
   [3]                  2     DDX11L1-202            havana
   [4]                  2     DDX11L1-202            havana
   [5]                  2     DDX11L1-202            havana
   [6]                  2     DDX11L1-201            havana
   [7]                  2     DDX11L1-201            havana
   [8]                  2     DDX11L1-201            havana
   [9]                  2     DDX11L1-201            havana
  [10]                  2     DDX11L1-201            havana
                       transcript_biotype         tag
                              <character> <character>
   [1]                               <NA>        <NA>
   [2]               processed_transcript       basic
   [3]               processed_transcript       basic
   [4]               processed_transcript       basic
   [5]               processed_transcript       basic
   [6] transcribed_unprocessed_pseudogene       basic
   [7] transcribed_unprocessed_pseudogene       basic
   [8] transcribed_unprocessed_pseudogene       basic
   [9] transcribed_unprocessed_pseudogene       basic
  [10] transcribed_unprocessed_pseudogene       basic
       transcript_support_level exon_number         exon_id exon_version
                    <character> <character>     <character>  <character>
   [1]                     <NA>        <NA>            <NA>         <NA>
   [2]                        1        <NA>            <NA>         <NA>
   [3]                        1           1 ENSE00002234944            1
   [4]                        1           2 ENSE00003582793            1
   [5]                        1           3 ENSE00002312635            1
   [6]                       NA        <NA>            <NA>         <NA>
   [7]                       NA           1 ENSE00001948541            1
   [8]                       NA           2 ENSE00001671638            2
   [9]                       NA           3 ENSE00001758273            2
  [10]                       NA           4 ENSE00001799933            2
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

• Another task: read the LEP (Leptin) transcript

# start your R session using 'R' command from the terminal
# 
rtracklayer::import(format="gtf", text=system("zcat Homo_sapiens.GRCh38.95.gtf.gz | grep -P 'gene_name \"LEP\"' ", intern=T))

Warning in readGFF(filepath, version = version, filter = filter): connection is not positioned at the start of the file, rewinding
  it

GRanges object with 12 ranges and 22 metadata columns:
       seqnames                 ranges strand |         source
          <Rle>              <IRanges>  <Rle> |       <factor>
   [1]        7 [128241284, 128257628]      + | ensembl_havana
   [2]        7 [128241284, 128257628]      + | ensembl_havana
   [3]        7 [128241284, 128241306]      + | ensembl_havana
   [4]        7 [128251991, 128252162]      + | ensembl_havana
   [5]        7 [128252019, 128252162]      + | ensembl_havana
   ...      ...                    ...    ... .            ...
   [8]        7 [128254404, 128254760]      + | ensembl_havana
   [9]        7 [128254761, 128254763]      + | ensembl_havana
  [10]        7 [128241284, 128241306]      + | ensembl_havana
  [11]        7 [128251991, 128252018]      + | ensembl_havana
  [12]        7 [128254764, 128257628]      + | ensembl_havana
                  type     score     phase         gene_id gene_version
              <factor> <numeric> <integer>     <character>  <character>
   [1]            gene      <NA>      <NA> ENSG00000174697            4
   [2]      transcript      <NA>      <NA> ENSG00000174697            4
   [3]            exon      <NA>      <NA> ENSG00000174697            4
   [4]            exon      <NA>      <NA> ENSG00000174697            4
   [5]             CDS      <NA>         0 ENSG00000174697            4
   ...             ...       ...       ...             ...          ...
   [8]             CDS      <NA>         0 ENSG00000174697            4
   [9]      stop_codon      <NA>         0 ENSG00000174697            4
  [10]  five_prime_utr      <NA>      <NA> ENSG00000174697            4
  [11]  five_prime_utr      <NA>      <NA> ENSG00000174697            4
  [12] three_prime_utr      <NA>      <NA> ENSG00000174697            4
         gene_name    gene_source   gene_biotype   transcript_id
       <character>    <character>    <character>     <character>
   [1]         LEP ensembl_havana protein_coding            <NA>
   [2]         LEP ensembl_havana protein_coding ENST00000308868
   [3]         LEP ensembl_havana protein_coding ENST00000308868
   [4]         LEP ensembl_havana protein_coding ENST00000308868
   [5]         LEP ensembl_havana protein_coding ENST00000308868
   ...         ...            ...            ...             ...
   [8]         LEP ensembl_havana protein_coding ENST00000308868
   [9]         LEP ensembl_havana protein_coding ENST00000308868
  [10]         LEP ensembl_havana protein_coding ENST00000308868
  [11]         LEP ensembl_havana protein_coding ENST00000308868
  [12]         LEP ensembl_havana protein_coding ENST00000308868
       transcript_version transcript_name transcript_source
              <character>     <character>       <character>
   [1]               <NA>            <NA>              <NA>
   [2]                  4         LEP-201    ensembl_havana
   [3]                  4         LEP-201    ensembl_havana
   [4]                  4         LEP-201    ensembl_havana
   [5]                  4         LEP-201    ensembl_havana
   ...                ...             ...               ...
   [8]                  4         LEP-201    ensembl_havana
   [9]                  4         LEP-201    ensembl_havana
  [10]                  4         LEP-201    ensembl_havana
  [11]                  4         LEP-201    ensembl_havana
  [12]                  4         LEP-201    ensembl_havana
       transcript_biotype         tag     ccds_id transcript_support_level
              <character> <character> <character>              <character>
   [1]               <NA>        <NA>        <NA>                     <NA>
   [2]     protein_coding       basic    CCDS5800                        1
   [3]     protein_coding       basic    CCDS5800                        1
   [4]     protein_coding       basic    CCDS5800                        1
   [5]     protein_coding       basic    CCDS5800                        1
   ...                ...         ...         ...                      ...
   [8]     protein_coding       basic    CCDS5800                        1
   [9]     protein_coding       basic    CCDS5800                        1
  [10]     protein_coding       basic    CCDS5800                        1
  [11]     protein_coding       basic    CCDS5800                        1
  [12]     protein_coding       basic    CCDS5800                        1
       exon_number         exon_id exon_version      protein_id
       <character>     <character>  <character>     <character>
   [1]        <NA>            <NA>         <NA>            <NA>
   [2]        <NA>            <NA>         <NA>            <NA>
   [3]           1 ENSE00001422382            2            <NA>
   [4]           2 ENSE00001303768            3            <NA>
   [5]           2            <NA>         <NA> ENSP00000312652
   ...         ...             ...          ...             ...
   [8]           3            <NA>         <NA> ENSP00000312652
   [9]           3            <NA>         <NA>            <NA>
  [10]        <NA>            <NA>         <NA>            <NA>
  [11]        <NA>            <NA>         <NA>            <NA>
  [12]        <NA>            <NA>         <NA>            <NA>
       protein_version
           <character>
   [1]            <NA>
   [2]            <NA>
   [3]            <NA>
   [4]            <NA>
   [5]               4
   ...             ...
   [8]               4
   [9]            <NA>
  [10]            <NA>
  [11]            <NA>
  [12]            <NA>
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

• You can, of course, read all of them them

# IT WILL TAKE LONG
foo<-rtracklayer::import("Homo_sapiens.GRCh38.95.gtf.gz")

SAM (BAM)

• https://en.wikipedia.org/wiki/SAM_(file_format)

• http://samtools.github.io/hts-specs/SAMv1.pdf

• You can check your bitwise FLAG interactively from here or here

• Let’s get the BAM file for SRR6481086 (you may need to install the sra toolkit)

sam-dump SRR6481086 | head -n 10

1   77  *   0   0   *   *   0   0   GGCCATGCAAGATTCCCATTCTTGCGACCCGGGTTCGTTTCCCGGGCGGCGCACCAGATCGGAAGAGCACACGNC AAAAAEEEA/EAEEEEE6EEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEE/EEEE6EEEE6AE6EEEEEEE#E RG:Z:TAATGCGC
1   141 *   0   0   *   *   0   0   GGTGCGCCGCCCNGGCCCCGCNNNCGNGTCGCCCGCCTNNNCNTCTTGCCTGNNCCGCTCGGCCGCGCGTCGTGT AAAAAEEEEEEE#EEA//EEE###EE#EEEE/EEEE/E###/#EAEEE//EE##AAEAE/EEA/EEE/EE/E/EA RG:Z:TAATGCGC
2   77  *   0   0   *   *   0   0   GAAGAGCACACGTCTGAACTCCAGTCACTAATGCGCATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAAAGGGNG AAAA/EE/E6EEEEEE66EEEE/EEE/EEE6EEEAE/A6EEEE6EEEEEE6EEAEE/EEEEEE6EEEEE//AE#E RG:Z:TAATGCGC
2   141 *   0   0   *   *   0   0   GGGGGGGGGGGGNGGGGGGGGNNNGGNGGGGGGGGGGGNNNGNGGGGGGGGGNNGGGGGGGGGGGGGGGGGGGGG AAAAAEEEEEEE#EEEEEEEE###EE#EEEEEEEEEEE###E#EEAEAEEEE##EEEEEEEEEEEEEAAAEEAAE RG:Z:TAATGCGC
3   77  *   0   0   *   *   0   0   GAAGAGCACACGTCTGAACTCCAGTCACTAATGCGCATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAAAGGGGG AAAA6EEAE6EEEEEE66EEEE/EEEEEE66EEEEE/EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE RG:Z:TAATGCGC
3   141 *   0   0   *   *   0   0   GGGGGGGGGGGGGGGGGGGGGNGNGGNGGGGGGGGGGGGGNGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG AAAAAEEEEEEEEEEEEEEEE#E#EE#EEEEEEEEEEEEE#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE RG:Z:TAATGCGC
4   77  *   0   0   *   *   0   0   GTAGAGCACACGTCTGAACTCCAGTCACTAATGCGCCTCTCGTATGCCGTCTTCTGCTTGAAAAAAAAAAGGGGG AAAA/AAEEEEEEEAE//EEAE/AEE6EEE/EEEEE/EEE/EE6EEEEEEEEEEAE/6EEEEE<EEEEE<6E6/E RG:Z:TAATGAGC
4   141 *   0   0   *   *   0   0   GGGGGGGGGGGGGGGGGGGGGGGNGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG AAAAAEEEEAEEEAEEEEEEEEE#EEEEEEEAEAAAEEEEEAAAAEEAAAAEAEAEAA//EAEAAE//AAEEEEA RG:Z:TAATGAGC
5   77  *   0   0   *   *   0   0   GAAGAGCACACGTCTGAACTCCAGTCACTAATGCGCATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAAAGGGGG AAAA/EE/E6EEEEEE//EEEEAEEE6EEAAEEEEE6EEEEEEEEEEEEEEEEEEE6EEEEEEEEEEEE6EEEEE RG:Z:TAATGCGC
5   141 *   0   0   *   *   0   0   GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG AAAAAEEEEEEEEAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE6EEEEEEEE RG:Z:TAATGCGC

• However, they turned out to be unaligned

# DO NOT RUN (takes a while)
sam-dump SRR6481086 | awk '{print $3}' | sort | uniq -c

• So, 1K reads of SRR6481086 were aligned to GRCh38 using HiSat2 (you are NOT expected to run this command)

# this is for your information only
# the path to the index should be changed (but I am not covering this)
hisat2 -x /home/ssg29/data/genome/Homo_sapiens/Ensembl/GRCh38/Sequence/HiSat2Index/genome --add-chrname --rg-id SRR6481086 --rg SRR6481086 -1 SRR6481086_pass_1.fastq.gz -2 SRR6481086_pass_2.fastq.gz | samtools view -bS - > SRR6481086.bam

• The BAM file is avaiable:

wget https://www.dropbox.com/s/6ouvliibmwhnbet/SRR6481086.bam

• Let’s see the header of the BAM file using samtools

# I assume samtools already installed 
# and the binary is avaialble from your $PATH environment
samtools view -H SRR6481086.bam

@HD VN:1.0  SO:unsorted
@SQ SN:chr1 LN:248956422
@SQ SN:chr2 LN:242193529
@SQ SN:chr3 LN:198295559
@SQ SN:chr4 LN:190214555
@SQ SN:chr5 LN:181538259
@SQ SN:chr6 LN:170805979
@SQ SN:chr7 LN:159345973
@SQ SN:chr8 LN:145138636
@SQ SN:chr9 LN:138394717
@SQ SN:chr10    LN:133797422
@SQ SN:chr11    LN:135086622
@SQ SN:chr12    LN:133275309
@SQ SN:chr13    LN:114364328
@SQ SN:chr14    LN:107043718
@SQ SN:chr15    LN:101991189
@SQ SN:chr16    LN:90338345
@SQ SN:chr17    LN:83257441
@SQ SN:chr18    LN:80373285
@SQ SN:chr19    LN:58617616
@SQ SN:chr20    LN:64444167
@SQ SN:chr21    LN:46709983
@SQ SN:chr22    LN:50818468
@SQ SN:chrX LN:156040895
@SQ SN:chrY LN:57227415
@SQ SN:chrMT    LN:16569
@RG ID:SRR6481086   SRR6481086
@PG ID:hisat2   PN:hisat2   VN:2.0.4    CL:"/usr/local/Cluster-Apps/hisat2/2.0.4/hisat2-align-s --wrapper basic-0 -x /home/ssg29/data/genome/Homo_sapiens/Ensembl/GRCh38/Sequence/HiSat2Index/genome --add-chrname --rg-id SRR6481086 --rg SRR6481086 -1 /tmp/89318.inpipe1 -2 /tmp/89318.inpipe2"

• Run samtools with ‘flagstat’

samtools flagstat SRR6481086.bam

2243 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
433 + 0 mapped (19.30%:-nan%)
2243 + 0 paired in sequencing
1242 + 0 read1
1001 + 0 read2
0 + 0 properly paired (0.00%:-nan%)
2 + 0 with itself and mate mapped
431 + 0 singletons (19.22%:-nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

• 2243 QC-passed reads?

samtools view -F 0x200 -c SRR6481086.bam

2243

• 0 QC-failed reads?

samtools view -f 0x200 -c SRR6481086.bam

0

• Why >1000*2 reads? (input was 2K)
• Hint: secondary alignment

samtools view -f 256 -c SRR6481086.bam

243

• For example:

samtools view SRR6481086.bam  | awk '$1=="SRR6481086.10"{print $0}'

SRR6481086.10   73  chr21   8259166 0   3S72M   =   8259166 0   GGCCTCCTCGTGGGGGGGCCGGGCCCCCCCTCCCCCGGCGCGACCGCTCTCCCACCCCTCCTCCCCGCGCCCCCG AAAAAEEEEEEEEE66//EEAE/EE6EEEE<EEE6E66EEE6/EE6E6E6EEE6EEEEAEEAEEEE/E6EEEEE/ AS:i:-11    ZS:i:-11    XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:22A8A40    YT:Z:UP RG:Z:SRR6481086 NH:i:3
SRR6481086.10   329 chr21   8442201 0   3S72M   =   8442201 0   GGCCTCCTCGTGGGGGGGCCGGGCCCCCCCTCCCCCGGCGCGACCGCTCTCCCACCCCTCCTCCCCGCGCCCCCG AAAAAEEEEEEEEE66//EEAE/EE6EEEE<EEE6E66EEE6/EE6E6E6EEE6EEEEAEEAEEEE/E6EEEEE/ AS:i:-11    ZS:i:-11    XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:22A8A40    YT:Z:UP RG:Z:SRR6481086 NH:i:3
SRR6481086.10   329 chr21   8214937 0   3S72M   =   8214937 0   GGCCTCCTCGTGGGGGGGCCGGGCCCCCCCTCCCCCGGCGCGACCGCTCTCCCACCCCTCCTCCCCGCGCCCCCG AAAAAEEEEEEEEE66//EEAE/EE6EEEE<EEE6E66EEE6/EE6E6E6EEE6EEEEAEEAEEEE/E6EEEEE/ AS:i:-11    ZS:i:-11    XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:22A8A40    YT:Z:UP RG:Z:SRR6481086 NH:i:3
SRR6481086.10   133 chr21   8259166 0   *   =   8259166 0   GGGGCGGGAGGGACGCCGCCGTCGCCGCCGCCCCCGAGCGCACGCTCCGCCGTCCCCCACAACGGGGGCCGCGCG AAA/AEEE/EAAA/EEEAEEAEEAEEEEAE/EEAAEAEEE///EAEAEAAEEAAAAAA/A//AA/A/AAAA/AEA YT:Z:UP RG:Z:SRR6481086

• Non-secondary alignment

samtools view -F 256 -c SRR6481086.bam

2000

• 443 mapped reads?

samtools view -f1 -F 0x4 -c SRR6481086.bam 

433

• 1242 read1 (from the 2243 paired in sequencing)?

samtools view -f 0x41 -c SRR6481086.bam 

1242

• 1001 read2 (from the 2243 paired in sequencing)?

samtools view -f 0x81 -c SRR6481086.bam 

1001

• 0 properly paired?
• Hint: both 0x1 and 0x2 bits set and 0x4 bit not set

samtools view -f 3 -F 0x4 -c SRR6481086.bam 

0

• 2 with itself and mate mapped
• Hint: exclude segment unmapped (0x4) AND next segment in the template unmapped (0x8)

samtools view -f 1 -F 12 -c SRR6481086.bam 

2

BED

• https://genome.ucsc.edu/FAQ/FAQformat#format1

• Let’s convert BAM format to BED format (using bedtools)

# I assume bedtools is exported to your $PATH
bedtools bamtobed -i SRR6481086.bam | head
# this is equivalent with the follwoing command
bamToBed -i SRR6481086.bam | head

chr21   8259165 8259237 SRR6481086.10/1 0   +
chr21   8442200 8442272 SRR6481086.10/1 0   +
chr21   8214936 8215008 SRR6481086.10/1 0   +
chr11   61967559    61967631    SRR6481086.11/1 60  -
chr2    32916247    32916322    SRR6481086.18/2 60  +
chrMT   2824    2896    SRR6481086.24/1 60  +
chr21   8442555 8442623 SRR6481086.27/1 0   +
chr21   8215292 8215360 SRR6481086.27/1 0   +
chr21   8398326 8398394 SRR6481086.27/1 0   +
chr21   8259521 8259589 SRR6481086.27/1 0   +
chr21   8259165 8259237 SRR6481086.10/1 0   +
chr21   8442200 8442272 SRR6481086.10/1 0   +
chr21   8214936 8215008 SRR6481086.10/1 0   +
chr11   61967559    61967631    SRR6481086.11/1 60  -
chr2    32916247    32916322    SRR6481086.18/2 60  +
chrMT   2824    2896    SRR6481086.24/1 60  +
chr21   8442555 8442623 SRR6481086.27/1 0   +
chr21   8215292 8215360 SRR6481086.27/1 0   +
chr21   8398326 8398394 SRR6481086.27/1 0   +
chr21   8259521 8259589 SRR6481086.27/1 0   +

• Count the mapped reads (in different ways)

bamToBed -i SRR6481086.bam | wc -l
samtools view -f 1 -F 4 SRR6481086.bam  | wc -l

433
433

• Compare the 0-based (e.g. BED) and 1-based (e.g. SAM) coordiates

samtools view -f 1 -F 4 SRR6481086.bam  | head | cut -f1-4 

SRR6481086.10   73  chr21   8259166
SRR6481086.10   329 chr21   8442201
SRR6481086.10   329 chr21   8214937
SRR6481086.11   89  chr11   61967560
SRR6481086.18   137 chr2    32916248
SRR6481086.24   73  chrMT   2825
SRR6481086.27   73  chr21   8442556
SRR6481086.27   329 chr21   8215293
SRR6481086.27   329 chr21   8398327
SRR6481086.27   329 chr21   8259522

• In theory, the aligned region can not go beyond > 75 for SRR6481086 (PE75), but there are:

bamToBed -i SRR6481086.bam | awk '$3-$2>75{print $0}'

chr7    32489330    32490374    SRR6481086.72/1 60  -
chr1    183627540   183630478   SRR6481086.78/1 60  -
chr19   49497979    49499519    SRR6481086.98/1 60  -
chr12   56713130    56713548    SRR6481086.173/1    60  -
chr20   3024592 3026730 SRR6481086.389/1    60  +
chr1    25229096    25232148    SRR6481086.493/1    60  -
chr16   81696643    81699732    SRR6481086.591/1    60  +
chr17   68524932    68525799    SRR6481086.846/1    60  +
chr22   39318550    39319623    SRR6481086.867/1    60  -

• Let’ look at the BAM file for those cases above

samtools view -f 1 -F 4 SRR6481086.bam  | awk '$1=="SRR6481086.72" || $1=="SRR6481086.78" || $1=="SRR6481086.98"{print $0}'

SRR6481086.72   89  chr7    32489331    60  14M975N55M6S    =   32489331    0   TTTGTCCACAAGCTCTAAGGGCAGCAGCTGCGACGGGTTGGTAGTAGCGTTAGCCGCCATGGCTACGCCCTCCGC EEEEAEEEEEEEEAEEEEEEEEEEEEEE6EAEEEEEE/6EE6EE6EEEEE/AEEE6EAE/EEE/A<////6AAAA AS:i:-6 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:69 YT:Z:UP RG:Z:SRR6481086 XS:A:-  NH:i:1
SRR6481086.78   89  chr1    183627541   60  54M2866N18M3S   =   183627541   0   TTTTCTTGCAGTCAAGACACGAACAATGGACCCTACTCCTCCAGCAGCAAGTGCCTTTTCATGCCATTGCAGGCC EEEE/EEA///E/E6A66E//666EEEA/E666E/6E66E666/66A666EE666EEEE666E//EEE//AAAAA AS:i:-3 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:72 YT:Z:UP RG:Z:SRR6481086 XS:A:-  NH:i:1
SRR6481086.98   89  chr19   49497980    60  67M1465N8M  =   49497980    0   TCCGCAAGTACAACCGCTTCGAGAAGCGCCACAAGAACATGTCTGTACACCTGTCCCCCTGCTTCAGGGACGTCC EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE/E/EEEEEEEEE6EEEEEEEE66EEE6EEAAAAA AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:75 YT:Z:UP RG:Z:SRR6481086 XS:A:+  NH:i:1

• You may see ‘N’ in CIGAR line, which mean ‘skipped region from the reference’
• This mean ‘intron-spanning’ reads in RNA-Seq (i.e. reads containing exon-intron junctions)
• See the GATK best practice with regard to this issue
• You can ‘split’ those N-containing CIGAR lines using bedtools

bamToBed -split -i SRR6481086.bam | awk '$3-$2>75{print $0}'

• Check out how the use of ‘-split’ option affects reporting overlap with the genome

Concluding Remarks

• There is more than one way to do it. (TMTOWTDI)
• You are not alone.
  • http://seqanswers.com/
  • https://www.biostars.org/
• It takes time.
• Good luck!