FastaAlternateReferenceMaker Consensus sequence header

onuralponuralp Barcelona, SpainMember
edited November 2013 in Ask the GATK team

Hi,

I am using FastaAlternateReferenceMaker to create consensus sequences as follows:

java -Xmx12g -jar /GenomeAnalysisTK-2.7-4-g6f46d11/GenomeAnalysisTK.jar -T FastaAlternateReferenceMaker -R /Reference/chromosome.1.fa -o /Output/consensus.1.fa --variant /VCF/chromosome1.vcf -L /Interval/chromosome.1.list

This command line produces a single fasta file including consensus sequences for the intervals provided in the list file.

Two questions:

(1) Presumably, sequence header for each sequence in the consensus fasta file corresponds to the contig name in the master reference file. Is there any way to modify the command line to print, say, the interval as the sequence header instead of the contig name?

Note: I can feed the program one interval at a time and name the output file accordingly, however, I'd like to stick to submitting my query per one chromosome at a time for the sake of saving up time.

(2) Number of sequences in the consensus fasta do not match the number of non-overlapping intervals in the input list. I would think that this is because some intervals are variant-free and therefore no alternate reference is reported for them. Do you confirm this is the case? I cannot easily check because the issue mentioned in (1).

Thank you.

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi there,

    1. It is not currently possible to create/name new contigs based on the intervals, no. The tool was not designed to do that.

    2. That should not be the case; in loci where no variant given, the reference base should still be emitted normally. How did you count the sequence lengths?

  • onuralponuralp Barcelona, SpainMember
    1. Ok, thank you.

    2. I have a list of genes with their chromosome start / end positions included in the interval list file.

    1:69090-70008 1:861120-879961 1:879582-894679 1:895966-901099 1:901876-910484 1:910578-917473 1:934341-935552 [...] 1:1017197-1051736 1:1072396-1079434 1:1108435-1114935 1:1109285-1133313 1:1138887-1142089 1:1146705-1149548

    I simply grepped the > character to count the number of sequences in the output fasta file.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Ah, I see. Keep in mind that if any of the intervals overlap (which includes intervals that don't share any positions but "touch", e.g. 1-9, 10-15) will be merged into a single interval.

  • onuralponuralp Barcelona, SpainMember

    Oh, that's indeed good to keep in mind - thank you.

Sign In or Register to comment.