Hi all! I am having problems and I hope I can get some help from you.

I will explain my situation: I'm trying to perform a PCA analysis to see how different several bam files are. I'm using the next pipeline:

 1. 
Getting the accession files. I am using the R library "SRAdb", so I am getting 4 files in .sra format.
 2. I use SRA-tools in order to convert the .sra file into .bam format
    with the following code:

     `sam-dump -r --min-mapq 25 $file | samtools view -bS > $file.bam`

 3. Sort `samtools sort $file -o $file_sorted`

 4. Index `samtools index $file_sorted $file_sorted.bai `

 5. Compute a matrix to generate the PCA plot

 ` multiBamSummary bins -b $files.bam -o my/out/path --smartLabels -bs 10000 -p 2`

At this point I'm getting the following error: 

    The file < myfile > does not have BAM or CRAM format. 

I haven't been able to trace the error, as any of the earlier steps reported any source of error. Any suggestions? (ideally I would like to skip the alignment step, I want to keep the file as original as possible)

 - sra-tools --version 2.9.1_1
 - samtools --version 1.9
 - deeptools --version 3.3.0

Thanks before hand!!


Let us use one of the example accession numbers above (SRR2089860). These are single-end reads.

Your options are:

Use `fastq-dump` to dump the reads out in fastq format (remove `-X 5` for full set)

    $ fastq-dump -X 5 SRR2089860
    Read 5 spots for SRR2089860
    Written 5 spots for SRR2089860

Use `sam-dump` to create fastq format files

    $ sam-dump --fastq SRR2089860 | head -16
    @HWI-D00473:169:HFK7WADXX:1:1101:1202:2011/1 unaligned
    NGAGTCTATACTCGTTACATTCGCGTAACTCATTGTTAATCGCGAAGTTGA
    +
    #1=DDDDFGHHGHJJJIJIIJGIIJJIJHICGIIIJJJIJGIJEHJIGIIG
    @HWI-D00473:169:HFK7WADXX:1:1101:1195:2074/1 unaligned
    CTCGAACTCCTCGTAGTGGCGATTGTCGGTGCTGCCCACCAGGTCCACTGT
    +
    CCCFFFFFHGHHHJIJIJJJJHIIGGGHIECEHFHGIEFIGGJGHJIIGIG
    @HWI-D00473:169:HFK7WADXX:1:1101:1230:2087/1 unaligned
    TGCCGGGAATTGTACAGTGCTCAGCTTTATAGGACATTTCCAAACAGTTAT
    +
    BBBFFFF8FHHHHJJJIJJJJIJGJJIJFGJIFGIIIJJJIGIEIIIIJGG
    @HWI-D00473:169:HFK7WADXX:1:1101:1222:2168/1 unaligned
    CCGAGACTTGCCTGCTCACCAGCGAAGAGGGCGAGGAGCGTTTGACGGCCG
    +
    @@CDDADDHFHHHIIIIIHGGE<GEGIEHIGIIDHGHGGIHHHEFFFCCCB

Use `sam-dump` to write SAM format files. This data appears to be unaligned (so --min-mapq should not affect anything, you can check).

    $ sam-dump -r SRR2089860 | head -4
    HWI-D00473:169:HFK7WADXX:1:1101:1202:2011       4       *       0       0       *       *       0       0       NGAGTCTATACTCGTTACATTCGCGTAACTCATTGTTAATCGCGAAGTTGA     #1=DDDDFGHHGHJJJIJIIJGIIJJIJHICGIIIJJJIJGIJEHJIGIIG
    HWI-D00473:169:HFK7WADXX:1:1101:1195:2074       4       *       0       0       *       *       0       0       CTCGAACTCCTCGTAGTGGCGATTGTCGGTGCTGCCCACCAGGTCCACTGT     CCCFFFFFHGHHHJIJIJJJJHIIGGGHIECEHFHGIEFIGGJGHJIIGIG
    HWI-D00473:169:HFK7WADXX:1:1101:1230:2087       4       *       0       0       *       *       0       0       TGCCGGGAATTGTACAGTGCTCAGCTTTATAGGACATTTCCAAACAGTTAT     BBBFFFF8FHHHHJJJIJJJJIJGJJIJFGJIFGIIIJJJIGIEIIIIJGG
    HWI-D00473:169:HFK7WADXX:1:1101:1222:2168       4       *       0       0       *       *       0       0       CCGAGACTTGCCTGCTCACCAGCGAAGAGGGCGAGGAGCGTTTGACGGCCG     @@CDDADDHFHHHIIIIIHGGE<GEGIEHIGIIDHGHGGIHHHEFFFCCCB

To do PCA analysis you will need to align fastq data to reference, count aligned reads to get an expression estimate. You could also use something like `salmon` to align to transcriptome to get counts.

