Skip to content

Instantly share code, notes, and snippets.

@jelber2
Created April 30, 2024 08:56
Show Gist options
  • Save jelber2/99d947ece717700bbe6f7cec9295381a to your computer and use it in GitHub Desktop.
Save jelber2/99d947ece717700bbe6f7cec9295381a to your computer and use it in GitHub Desktop.
Converts quality scores in BAM alignments to a single ASCII Phred Score you desire
# note: outputs to STDOUT a SAM file without a header
# note2: remove auxillary tags and information
# note3: ignores RNEXT and PNEXT from BAM file (puts * and 0 respectivelv)
# tested with julialang v1.10.2 and XAM v0.4.0
# $ julia changeBAMquality.jl input.bam ? > test.sam.noheader
#
# to add a header from the original BAM file that had its qualities changed
# and add on a PG line
# $ ASCII_CHARACTER="?"
# $ INPUT=input.bam
# $ INPUT2=test.sam.noheader
# this is the command to make the BAM file
# $ cat <(samtools view -H $INPUT) <(echo -e "@PG\tID:changeBAMquality\tPN:changeBAMquality.jl\tVN:0.01\tCL:julia changeBAMquality.jl $INPUT $ASCII_CHARACTER") $INPUT2 |samtools view -Sb > test.bam
#
using XAM
using Base.Iterators: partition
using ArgParse
function parse_commandline()
s = ArgParseSettings()
@add_arg_table s begin
"input_file"
help = "Path to the input BAM file"
required = true
"ASCII_character"
help = "ASCII character quality value to change, for examle ? for Q30"
required = true
arg_type = String
end
return parse_args(s)
end
function change_quality(input_file::String, ASCII_character::String)
reader = open(BAM.Reader, input_file)
record = BAM.Record()
while !eof(reader)
empty!(record)
read!(reader, record)
if BAM.ismapped(record)
if BAM.isprimaryalignment(record)
println(BAM.tempname(record), "\t", BAM.flags(record), "\t", BAM.refname(record), "\t", BAM.position(record), "\t", BAM.mappingquality(record), "\t", BAM.cigar(record), "\t*", "\t0\t", BAM.templength(record), "\t", BAM.sequence(record), "\t", repeat(ASCII_character, BAM.seqlength(record)))
end
end
end
close(reader)
end
function main()
parsed_args = parse_commandline()
input_file = parsed_args["input_file"]
ASCII_character = parsed_args["ASCII_character"]
if length(ASCII_character) != 1
println("Error: Input must be a single character.")
return
end
if !isascii(ASCII_character[1])
println("Error: Input must be an ASCII character.")
return
end
change_quality(input_file, ASCII_character)
end
main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment