1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
| using BioSequences
using FASTX
function cal(sequence::LongDNASeq)
nGC = count(x -> ((x == DNA_G) || (x == DNA_C)), sequence)
nAT = count(x -> ((x == DNA_A) || (x == DNA_T)), sequence)
(nGC, nAT)
end
function getPosGCs(IO::IOStream, CHRseq::BioSequence, totalLen::Int64, Chr::SubString, POS::Int64, WINDOW::Vector{Int64}, THRESH::Int64, idx::SubString)
print(IO, idx, "\t", Chr, "\t", POS, "\t")
for window in WINDOW
window = window % 2 == 0 ? window + 1 : window
startPOS = POS - Int64(floor(window / 2))
tailPOS = POS + Int64(floor(window / 2))
tailPOS = tailPOS > totalLen ? totalLen : tailPOS
startPOS = startPOS <= 0 ? 1 : startPOS
gc, at = cal(CHRseq[startPOS:tailPOS])
if gc + at > THRESH
print(IO, round(gc / (gc + at); digits=6), "\t")
else
print(IO, "NA", "\t")
end
end
print(IO, "\n")
end
function main(snpLoci::String, fasta::String)
open("GC_correct.file.txt", "w") do io
println("Loding reference genome!")
Genome = open(FASTA.Reader, fasta, index=string(fasta, ".fai"))
WINDOWS = Int64[25, 50, 100, 200, 500, 1e3, 2e3, 5e3, 1e4, 2e4, 5e4, 1e5, 2e5, 5e5, 1e6]
println(io, "\tChr\tPosition\t25bp\t50bp\t100bp\t200bp\t500bp\t1kb\t2kb\t5kb\t10kb\t20kb\t50kb\t100kb\t200kb\t500kb\t1Mb")
chrMarker = ""
seq = dna"NNNNN"
totalLen = 0
for line in eachline(snpLoci)
if startswith(line, "\t")
continue
else
idx, chr, pos = split(line, "\t")
pos = parse(Int64, pos)
if chr != chrMarker
seq = FASTA.sequence(Genome[chr])
totalLen = length(seq)
println("processing chromosome ", chr, "...")
end
chrMarker = chr
getPosGCs(io, seq, totalLen, chr, pos, WINDOWS, 20, idx)
end
end
end
end
main(ARGS[1], ARGS[2])
|
