...
The read files were downloaded from the ENA SRA study.
So that we can treat all the data as single-ended for simplicity, we concatenated two separate FASTQ (paired-end) files for sample SRR030252 using this command
Code Block |
---|
cat SRR030252_1.fastq SRR030252_2.fastq > SRR030252.fastq
|
The reference genome file was downloaded from the NCBI Genomes page.
We renamed the FASTA sequence from "gi|254160123|ref|NC_012967.1|" to "NC_012967" by changing the first line of the NC_012967.1.fasta sequence using a text editor.
Map Reads
Choose an appropriate program and map the reads. As for other variant callers, convert the mapped reads to BAM format, then sort and index the BAM file.
...
Code Block | ||
---|---|---|
| ||
login1$ R ... > regions <- data.frame(chr="gi|254160123|ref|NC_012967.1|", start = 1, stop=100000) > result = deepSNV(test = "SRR032376.sorted.bam", control = "SRR030252.sorted.bam", regions=regions) > sig_result <- summary(result, sig.level=0.05, adjust.method="BH") > pdf("SRR032374.pdf") > plot(mix) > dev.off() > write.csv(sig_result, "SRR032374") |
...