I am developing a pipeline (with GATK v3.7) for data from whole-exome sequencing from Illumina NextSeq 500. I have FASTQ files, and I want to create a very small dummy dataset from this, so that I can easily test tweaks to my pipeline, without having to wait for hours every time I tweak something.
Can I simply use, say, the first 1000 reads from a file, and expect that to work as a test dataset?
2 answers
Can I simply use, say, the first 1000 reads from a file
Don't do that especially since the first reads in an illumina dataset may not be of good quality. You can use reformat.sh from BBMap suite to sample reads (options below) in various ways (by number of reads, bases, a certain fraction (which may be best)).
reads=-1 Set to a positive number to only process this many INPUT reads (or pairs), then quit.
skipreads=-1 Skip (discard) this many INPUT reads before processing the rest.
samplerate=1 Randomly output only this fraction of reads; 1 means sampling is disabled.
sampleseed=-1 Set to a positive number to use that prng seed for sampling (allowing deterministic sampling).
samplereadstarget=0 (srt) Exact number of OUTPUT reads (or pairs) desired.
samplebasestarget=0 (sbt) Exact number of OUTPUT bases desired.
BTW: This would not be a dummy dataset but a representative one.
For testing a variant-calling pipeline, you may want sufficient coverage to actually call a variant. To do so I would suggest selecting some 100kbp chunk of the genome (which should be enough for perhaps 100 variants, including indels and substitutions) and pulling out all the reads that map to it. Then use those reads as your representative dataset (and you can use the 100kbp genome chunk as your representative reference); it will have sufficient coverage to actually call variants normally, but will still be a very small file.
I don't see how I might do this
You will start with your aligned bam file and then use (Extract Reads From A Bam File That Fall Within A Given Region ) to extract reads from a particular region you are interested in.
Ah, I see. Using: samtools view -b input.bam "chr10:1-10000000" > output.bam, for example. This is a good idea.
It presents some challenges in my particular case, however. I need separate FASTQs for paired-end reads and four lanes to test my pipeline, and while I can extract the four lanes using the read group information, I can't do this for the paired-end reads.
What I typically do is open a fasta file in a text editor and manually select a bunch of text and paste it into a new fasta file (so that I can make sure I'm not selecting 100kbp of Ns or something). Then:
bbmap.sh in1=r1.fq in2=r2.fq outm1=mapped1.fq outm2=mapped2.fq ref=100kbp.fa
Log in to answer this question.
-
Illumina simulated reads - tools or SRA projects
written by aziznasr1920 •
-
Does GATK provide datasets to test their pipelines? (Learning purpose)
written by yussab •
-
Analysis of exome samples build from different exome kits
written by Nicolas Rosewick
-
Sample Mouse RNA-Seq Dataset
written by isu2017 •
-
Find illumina adaptors
written by sabaghianamir70 •
-
Test set for WGS analysis pipeline
written by gwotto •
-
Is it worth realigning 1000 Genomes?
written by Jeremy Leipzig
-
As an illumina customer, do you use basespace, alternative cloud solution or pipelines built in-hou…
written by mani824 •
-
Gatk Snp Calling Test Set
written by Tryingtogetthere •
-
Analysing Multiple Exome Datasets
written by User 9126 •
you can do that or you can random sample your fastq.
For first 1000 sequences:
For random sampling 1000 reads from your fastq:
Download seqkit from here: http://bioinf.shenwei.me/seqkit/download/
small typo in the second command