2.9 years ago by
This command in the User Guide uses the program
samtools faidx to fetch the guide sequence from the reference genome. The reference genome should be in
.fasta format. In your command above, you are using a
.bam file of mapped reads rather than the reference genome. The command is also missing the call to
samtools. So, if your reference genome is called
genome.fasta, the command should look like this:
reference=system(sprintf("samtools faidx genome.fasta %s:%s-%s", seqnames(gdl), start(gdl), end(gdl)), intern = TRUE)
For this to work, you will need to have
samtools installed. You can use the R function
file.exists() to check if your filenames are correct.
If you have mapped your reads to a single amplicon sequence rather than a full genome, an alternative to using
samtools is to read in the amplicon sequence using
Biostrings then select the reference sequence using
amplicon <- Biostrings::readDNAStringSet("genome.fasta")
reference <- substr(amplicon, guide_start, guide_end)
Or if you already know the sequence of the region you want to use, you could create the reference sequence directly, e.g.
reference <- Biostrings::DNAString("AACCTTGG")
Hope this helps, and happy new year!