2
$\begingroup$

What's the easiest way to run Picard GCBias ignoring certain chromosomes/sequences in the reference? Looking at the CollectGcBiasMetrics, it seems there isn't a bed file option that can be passed.

java -jar ~/picard-tools-1.134/picard.jar CollectGcBiasMetrics --help
USAGE: CollectGcBiasMetrics [options]

Documentation: http://broadinstitute.github.io/picard/command-line-overview.html#CollectGcBiasMetrics

Tool to collect information about GC bias in the reads in a given BAM file. Computes the number of windows (of size specified by WINDOW_SIZE) in the genome at each GC% and counts the number of read starts in each GC bin.  What is output and plotted is the "normalized coverage" in each bin - i.e. the number of reads per window normalized to the average number of reads per window across the whole genome..

Version: 1.134(a7a08c474e4d99346eec7a9956a8fe71943b5d80_1434033355)


Options:

--help
-h                            Displays options specific to this tool.

--stdhelp
-H                            Displays options specific to this tool AND options common to all Picard command line 
                              tools.

--version                     Displays program version.

CHART_OUTPUT=File
CHART=File                    The PDF file to render the chart to.  Required. 

SUMMARY_OUTPUT=File
S=File                        The text file to write summary metrics to.  Default value: null. 

WINDOW_SIZE=Integer           The size of windows on the genome that are used to bin reads.  Default value: 100. This 
                              option can be set to 'null' to clear the default value. 

MINIMUM_GENOME_FRACTION=DoubleFor summary metrics, exclude GC windows that include less than this fraction of the 
                              genome.  Default value: 1.0E-5. This option can be set to 'null' to clear the default 
                              value. 

IS_BISULFITE_SEQUENCED=Boolean
BS=Boolean                    Whether the SAM or BAM file consists of bisulfite sequenced reads.  Default value: false. 
                              This option can be set to 'null' to clear the default value. Possible values: {true, 
                              false} 

METRIC_ACCUMULATION_LEVEL=MetricAccumulationLevel
LEVEL=MetricAccumulationLevel The level(s) at which to accumulate metrics.  Default value: [ALL_READS]. This option can 
                              be set to 'null' to clear the default value. Possible values: {ALL_READS, SAMPLE, 
                              LIBRARY, READ_GROUP} This option may be specified 0 or more times. This option can be set 
                              to 'null' to clear the default list. 

INPUT=File
I=File                        Input SAM or BAM file.  Required. 

OUTPUT=File
O=File                        File to write the output to.  Required. 

ASSUME_SORTED=Boolean
AS=Boolean                    If true (default), then the sort order in the header file will be ignored.  Default 
                              value: true. This option can be set to 'null' to clear the default value. Possible 
                              values: {true, false} 

STOP_AFTER=Long               Stop after processing N reads, mainly for debugging.  Default value: 0. This option can 
                              be set to 'null' to clear the default value. 

Given this, and wanting to calculate GCBias the using Picard (not other software), what's the best option?

$\endgroup$

1 Answer 1

1
$\begingroup$

The proper solution is to use a different tool (some of this you could do with computeGCBias in deepTools). But since you don't want to do that, you'll have to manually remove unwanted chromosomes from the BAM file(s):

$ cat foo.awk
BEGIN{split(excludes, excludeList, " ")}
{exclude=0
if($1 == "@SQ") {
    for(ex in excludeList) {
        if(ex == substr($2, 4)) {
            exclude=1
            break
        }
    }
} else { \
    for(ex in excludeList) {
        if(ex == $3) {
            exclude=1
            break
        }
    }
}
if(!exclude) print}

$ excludes="chr1 chr2"
$ samtools view -h alignments.bam | awk -v excludes="$excludes" -f foo.awk | samtools -bo alignments.filtered.bam

The awk script is stripping unwanted chromosomes from the header and the body. You can then run the result through picard.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.