Page Comparison

#look at the set of barcodes included for our samples
lesscat barcodes_Lib1.tsv 


#check the reads
grep CATG Lib01_R1.fastq | head


#check where these are in our reads
grep "^TCGAT" Lib01_R1.fastq | wc -l
grep "^CGATC" Lib01_R1.fastq | wc -l
grep "^AACCA" Lib01_R1.fastq | wc -l
grep "^GCATG" Lib01_R1.fastq | wc -l

#make a directory to put the resulting sample fastq files into
mkdir sample_fastqs
 
#look at the documentation for process_radtags
./process_radtags -h
  
#execute the command to process the rad data
./process_radtags -i 'fastq' -1 Lib01_R1.fastq -2 Lib01_R2.fastq -o ./sample_fastqs/ -b barcodes_Lib1.tsv --inline_index -e 'nlaIII' -r --disable_rad_check

#check out our reference
less stickleback_chrom3.fasta 


#how many sequences are there?
grep "^>" stickleback_chrom3.fasta | wc -l


#index the reference
module load bwa
bwa index stickleback_chrom3.fasta


#run mapping
bwa mem stickleback_chrom3.fasta sample_fastqs/sample_AACCA-AGCGAC.1.fq sample_fastqs/sample_AACCA-AGCGAC.2.fq > sampleA.sam
bwa mem stickleback_chrom3.fasta sample_fastqs/sample_CGATC-AGCGAC.1.fq sample_fastqs/sample_CGATC-AGCGAC.2.fq > sampleB.sam
bwa mem stickleback_chrom3.fasta sample_fastqs/sample_GCATG-AGCGAC.1.fq sample_fastqs/sample_GCATG-AGCGAC.2.fq > sampleC.sam
bwa mem stickleback_chrom3.fasta sample_fastqs/sample_TCGAT-AGCGAC.1.fq sample_fastqs/sample_TCGAT-AGCGAC.2.fq > sampleD.sam


#assess mapping efficiency
module load samtools
samtools flagstat sampleA.sam > flagStatRes_sampleA.txt
samtools flagstat sampleB.sam > flagStatRes_sampleB.txt
samtools flagstat sampleC.sam > flagStatRes_sampleC.txt
samtools flagstat sampleD.sam > flagStatRes_sampleD.txt


#these results are not spectacular.
#80% mapping efficiency is OK, but the low rate of properly paired reads is disapointingdissapointing.


#convert the sam files to bam files
samtools sort sampleA.sam -O SAMBAM -o sampleA.bam
samtools sort sampleB.sam -O SAMBAM -o sampleB.bam
samtools sort sampleC.sam -O SAMBAM -o sampleC.bam
samtools sort sampleD.sam -O SAMBAM -o sampleD.bam

#index the reference for samtools
module load samtools
samtools faidx stickleback_chrom3.fasta
 
#make a list of the bam files
ls *.bam > my_bamfiles.txt
 
#run mpileup
samtools mpileup -f stickleback_chrom3.fasta -t DP,AD,ADF,ADR,SP -u -b my_bamfiles.txt > mpileup_results.bcf


#now call genotypes from the mpileup results
bcftools call -vmO v -o raw_calls.vcf mpileup_results.bcf