Once variant calling is completed, this will create a variant call format (VCF) file. NB: Variant calling takes about 5-10 minutes! Be patient. In your projects you will be dealing with fission yeast which is haploid. At these variant positions or sites, Freebayes will also genotype each of the three samples, i.e. at each of these positions, it will determine whether each sample is homozygous for the reference allele, heterozygous or homozygous for the alternate allele): freebayes -f TriTrypDB-39_LdonovaniBPK282A1_Genome.fasta -ploidy 2 -bam-list bams.list > & Here you are redirecting the output of the ls (list) command to a file called bams.listįreebayes will identify variants (SNPs and indels) that are different from the reference genome across your three samples. So you can tell FreeBayes which samples to use for variant calling, make a file containing a list of the relevant bam files: ls *.rmdup.bam > bams.list See this description to understand the fai file format. It contains information about the number and size of the chromosomes (strictly speaking the scaffolds making up the reference genome). Now index, pull out chromosome 1 reads, and mark PCR duplicates for samples 2 and 3.įirst ‘index’ the fasta format file of the genome with samtools samtools faidx TriTrypDB-39_LdonovaniBPK282A1_Genome.fasta Remember, to come out of less and return to the command line, press the q key. The 8th line of this file contains information about the percentage of duplication (second number from the end). Use less/more to have a look at the metrics file. You can do this using Picard Tools MarkDuplicates like so: java -XX:ParallelGCThreads=2 -Xmx2g -jar picard.jar MarkDuplicates \ The next step is to remove PCR duplicates. bai samtools index .bamĮxtract only sequence reads that have aligned to chromosome 1: samtools view -b .bam chr1 > .chr1.bam &Ĭheck that you have successfully created the chromosome 1 file, and find out the size of the file. This will create an index file for each BAM file ending in. This will allow you to complete the variant calling during the workshop.įirst index the three BAM files that you have just copied over. This avoids having to copy the picard program into your directory.įor this workshop you will only call variants on chromosome 1, which is 280kb long (the whole genome of Leishmania is about 30Mb long). If you list the files in your directory, you will see that you now have a link file called picard.jar that points to the location of the executable picard.jar file. ln -sf /opt/yarcc/data/biology/mbiol_sequence_analysis_2019/picard.jar. Make a symbolic link to the java executable picard.jar file. Picard can do a lot of other useful things too, see here. You will use the MarkDuplicates tool from Picard Tools, to mark PCR duplicates so that they are not used in the subsequent steps. It is important to remove/mark PCR duplicate sequences as they are non-independent. Load the modules (i.e. the programs that you will use) source /opt/yarcc/data/biology/mbiol_sequence_analysis_2019/biolprogramsĪs in the previous workshop, open a second connection to your designated server in which you can run top to keep an eye on tasks you are running: top -u userid Once you are logged into your designated server, use the cd command to navigate to the analysis directory you made last time. To do this, first login in to biollogin, and then connect to biolnode0: ssh biolnode0 If your userid is in the alphabetic range ag1314 to ja1082, you can run the analysis for the workshops (and your projects) on the server called .uk: ssh your username is in the alphabetic range jdh562 to xc1051, you will have you run your analysis for the workshops (and your projects) on a different server called .uk. Log on to your designated server as you did in Workshop 2 using putty (from a PC) or a terminal window (from a Mac).