Last active
May 2, 2018 23:45
-
-
Save jodyphelan/e721112d91dcaaa67cd79217ef2af24f to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| git clone https://github.com/lh3/seqtk.git && cd seqtk && make && cd .. | |
| wget https://github.com/samtools/bcftools/releases/download/1.8/bcftools-1.8.tar.bz2 && tar -xvf bcftools-1.8.tar.bz2 && cd bcftools-1.8 && make && cd .. | |
| wget https://github.com/samtools/samtools/releases/download/1.8/samtools-1.8.tar.bz2 && tar -xvf samtools-1.8.tar.bz2 && cd samtools-1.8/ && make && cd .. | |
| git clone https://github.com/lh3/minimap2.git && cd minimap2 && make && cd .. | |
| wget https://zlib.net/pigz/pigz-2.4.tar.gz && tar -xvf pigz-2.4.tar.gz && cd pigz-2.4/ && make && cd .. | |
| wget ftp://ftp.ensemblgenomes.org/pub/release-32/bacteria//fasta/bacteria_0_collection/mycobacterium_tuberculosis_h37rv/dna/Mycobacterium_tuberculosis_h37rv.ASM19595v2.dna.toplevel.fa.gz | |
| gunzip Mycobacterium_tuberculosis_h37rv.ASM19595v2.dna.toplevel.fa.gz | |
| mv Mycobacterium_tuberculosis_h37rv.ASM19595v2.dna.toplevel.fa H37Rv.fa | |
| wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR611/008/SRR6117388/SRR6117388.fastq.gz | |
| ./seqtk/seqtk sample -s 123 SRR6117388.fastq.gz 0.2 | ./pigz-2.4/pigz -c > SRR6117388.sampled.fastq.gz | |
| export THREADS=40 | |
| ./minimap2/minimap2 -t $THREADS -R '@RG\tID:test\tSM:test\tPL:minION' -a -x map-ont H37Rv.fa SRR6117388.sampled.fastq.gz | ./samtools-1.8/samtools view -@ $THREADS -b - | ./samtools-1.8/samtools sort -@ $THREADS -o test.bam - | |
| ./samtools-1.8/samtools index test.bam | |
| # Call variants in region and look at variant at position 3329088 | |
| ./bcftools-1.8/bcftools mpileup -f H37Rv.fa test.bam -r Chromosome:3320088-3330088 -B -Q8 -a AD | ./bcftools-1.8/bcftools call -mv | grep 3329088 | |
| # Lookin at the AD tag it looks like there are 3 reads supporting the ALT | |
| # Look at the mpileup from samtools using the same threshold | |
| ./samtools-1.8/samtools mpileup -f H37Rv.fa test.bam -r Chromosome:3328088-3330088 -B -Q8 | grep 3329088 | |
| # The large variant does not appear | |
| # Look at the mpileup from samtools using more relaxed Q threshold | |
| ./samtools-1.8/samtools mpileup -f H37Rv.fa test.bam -r Chromosome:3328088-3330088 -B -Q0 | grep 3329088 | |
| # Now the variant appears |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment