Graphing GC Bias in Sequence Capture

Question

[Jon17]

Anyone have recommendations for making a chart like this for sequence capture?

Picard CollectGcBiasMetrics only works for whole genome analysis.

Picard hsMetrics tends do do funny things with window sizes.

I've had some success with a bedtools & custom scripting but can't reproduce other people's material. I'm close, but I don't know if they are right or if I am right.

Also, anyone got any favorite illumina read simulators that creates reads given a bed file and produces a known gc bias? I'm trying to make control data.

Answer 1

[trausch]

Alfred calculates the reference GC, target region GC and sample GC content. Here is an example, just select "DNA-Seq, Exome" from the reference data set collection and then jump to "GC Content". Command should be:

alfred qc -b targets.bed -r genome.fa -o out.stats.tsv.gz -j out.stats.json.gz exome.bam

[Jon17]

Thank you!!!!!!!!!!!!!!!!!

Answer 2

[Chao Tang]

May be the function letterFrequency in package Biostrings is a good choice for GC content calculating.

gs_seq <- BSgenome::getSeq(Hsapiens, gr) letterFrequency(gs_seq, "GC", as.prob=TRUE)