Easiest way to get number of reads per contig?

[kmkocot]

Hi all,

I'm trying to map 2 X 100 PE Illumina reads (transcriptome) back to a small number of contigs. I want to get the number of reads that align to each contig as quickly as possible. Here's the commands I've tried so far but it isn't working.

SCAFFOLD_INPUT_FILE=test_selected_sequences.fas
S1_FASTQ_INPUT_FILE=D13F4ACXX_s6_1_GSLv2-7_32_SL18114_test.fastq.gz
S2_FASTQ_INPUT_FILE=D13F4ACXX_s6_2_GSLv2-7_32_SL18114_test.fastq.gz
REF_MAP_BASENAME=test

BWA=`which bwa`
SAMTOOLS=`which samtools`
BCFTOOLS=`which bcftools`
VCFUTILS=`which vcfutils.pl`
BT2=`which bowtie2`
BT_BUILD=`which bowtie2-build`
TDSTAMP=`date | awk 'OFS="_" {print ($2,$3,$6,$4)}'`

${BT_BUILD} ${SCAFFOLD_INPUT_FILE} ${REF_MAP_BASENAME}
${BT2} -x ${REF_MAP_BASENAME} -1 ${S1_FASTQ_INPUT_FILE} -2 ${S2_FASTQ_INPUT_FILE} --qc-filter --no-unal -S ${REF_MAP_BASENAME}.bowtie2.mapped.sam
$SAMTOOLS view -bT ${SCAFFOLD_INPUT_FILE} ${REF_MAP_BASENAME}.bowtie2.mapped.sam > ${REF_MAP_BASENAME}.bowtie2.mapped.bam
$SAMTOOLS sort ${REF_MAP_BASENAME}.bowtie2.mapped.bam ${REF_MAP_BASENAME}.bowtie2.mapped.sorted
$SAMTOOLS idxstats ${REF_MAP_BASENAME}.bowtie2.mapped.sorted.bam

Here's the errors I get:

[samopen] SAM header is present: 28 sequences.
[bam_index_load] fail to load BAM index.
[bam_idxstats] fail to load the index.

The problem seems to be once the samtools apps start running. I'm using samtools version 0.1.19-44428cd. Any ideas on what I'm doing wrong? Is there an easier way?

Thanks!
Kevin

[sdriscoll]

revise:

${BT_BUILD} ${SCAFFOLD_INPUT_FILE} ${REF_MAP_BASENAME}
${BT2} -x ${REF_MAP_BASENAME} -1 ${S1_FASTQ_INPUT_FILE} -2 ${S2_FASTQ_INPUT_FILE} --qc-filter --no-unal -S ${REF_MAP_BASENAME}.bowtie2.mapped.sam
$SAMTOOLS view -bT ${SCAFFOLD_INPUT_FILE} ${REF_MAP_BASENAME}.bowtie2.mapped.sam > ${REF_MAP_BASENAME}.bowtie2.mapped.bam
$SAMTOOLS sort ${REF_MAP_BASENAME}.bowtie2.mapped.bam ${REF_MAP_BASENAME}.bowtie2.mapped.sorted
$SAMTOOLS idxstats ${REF_MAP_BASENAME}.bowtie2.mapped.sorted.bam

to

${BT_BUILD} ${SCAFFOLD_INPUT_FILE} ${REF_MAP_BASENAME}
${BT2} -x ${REF_MAP_BASENAME} -1 ${S1_FASTQ_INPUT_FILE} -2 ${S2_FASTQ_INPUT_FILE} --qc-filter --no-unal -S ${REF_MAP_BASENAME}.bowtie2.mapped.sam
$SAMTOOLS view -bT ${SCAFFOLD_INPUT_FILE} ${REF_MAP_BASENAME}.bowtie2.mapped.sam > ${REF_MAP_BASENAME}.bowtie2.mapped.bam
$SAMTOOLS sort ${REF_MAP_BASENAME}.bowtie2.mapped.bam ${REF_MAP_BASENAME}.bowtie2.mapped.sorted
$SAMTOOLS index ${REF_MAP_BASENAME}.bowtie2.mapped.sorted.bam
$SAMTOOLS idxstats ${REF_MAP_BASENAME}.bowtie2.mapped.sorted.bam

I added in the call to "samtools index" which is why the samtools idxstats call is failing.

[kmkocot]

Thanks! I'm testing it now.