Advice on tools sought - exporting consensus sequences with coverage & SNPs?

[TabeaK]

Hi all!

I need some tools advice.
So; I have mapped reads in BAM format.
I have called SNPs using samtools/bcftools.

Are there any good tools for exporting consensus sequences from; say a vcf file; in; say fasta format; that maintain infos like coverage?

Hope I am making sense here. Feel free to point me to appropriate threads on seqanswers; my crappy searching has not turned up anything!

[TabeaK]

Okay: after a lot of reading around here on seqanswers, I think I have managed to some of this with samtools - always helps to read the manual carefully.

Only thing where I am stuck is the consensus. My current read library has okay-ish, but not great, coverage which means there are gaps in my alignment to the reference. These gaps might be there because of extensive deletion or something: but lets not even worry about that now...

So, if I call a consensus from the vcf file; how are the regions having no coverage treated? Is the sequence just taken from the ref sequence? Can I track/visualise that somehow? Just so I know which bit in the consensus came from the actual bam (i.e. has experimental evidence from my sequencing runs) and which bit came from the reference and has no experimental evidence?

[TabeaK]

Some progress.

Apologies if it seems like I am talking to myself; but maybe at some point someone will search for the same issue and this may help...

So; after using the vcfutils and f2q conversion in samtools; I have ended up with a consensus fastq file.

The file contains: Upper case letters (good quality?); lower case letters (low quality?) and - most interestingly - a bunch of Ns (lower case) where my reads do not cover the reference.

After converting this to fasta (seqtk) - which obviously looses the upper/lower case info - I have something I can work with reasonably well for now!

2 Questions remain for me:

1. Is there a description of what vcfutils actually does in detail somewhere? So, when does it call upper/lower case?
2. Is there any way to force it to use the reference sequence when it does not find coverage? So; it prints the ref sequence instead of an N?

Maybe I need another tool for that.

[swbarnes2]

Some progress.

Apologies if it seems like I am talking to myself; but maybe at some point someone will search for the same issue and this may help...

So; after using the vcfutils and f2q conversion in samtools; I have ended up with a consensus fastq file.

The file contains: Upper case letters (good quality?); lower case letters (low quality?) and - most interestingly - a bunch of Ns (lower case) where my reads do not cover the reference.

After converting this to fasta (seqtk) - which obviously looses the upper/lower case info - I have something I can work with reasonably well for now!

2 Questions remain for me:

1. Is there a description of what vcfutils actually does in detail somewhere? So, when does it call upper/lower case?
2. Is there any way to force it to use the reference sequence when it does not find coverage? So; it prints the ref sequence instead of an N?

Maybe I need another tool for that.

I have an older version of samtools but vcfutils is in perl, so it's human readable. If you are using vcf2fastq, the key lines are these:

my ($b, $q);
	  $q = $1 if ($t[7] =~ /FQ=(-?[\d\.]+)/);
	  if ($q < 0) {
		$_ = ($t[7] =~ /AF1=([\d\.]+)/)? $1 : 0;
		[B]$b = ($_ < .5 || $alt eq '.')? $ref : $alt;[/B]
		$q = -$q;
	  } else {
		[B]$b = $het{"$ref$alt"};
		$b ||= 'N';[/B]
	  }
	  $b = lc($b);
	  $b = uc($b) if (($t[7] =~ /MQ=(\d+)/ && $1 >= $_Q) && ($t[7] =~ /DP=(\d+)/ && $1 >= $_d && $1 <= $_D));

$b is depending on the vcf entry, either the reference letter, or the alt letter, or a letter representing both letters, if the vcf thinks its a mixed SNP, or an N, if there's nothing at all for that letter. Those lines are the first 3 bolded lines.

The last line determines if $b is uppercase or lowercase, which basically happens if the mapping quality is above the minimum you specified, and the coverage is between the min and max you specified.

The other place in the program where it will make your sequence lowercase is here:

for my $g (@$gaps) {
	my $beg = $g->[0] > $l? $g->[0] - $l : 0;
	my $end = $g->[0] + $g->[1] + $l;
	$end = length($$seq) if ($end > length($$seq));
	substr($$seq, $beg, $end - $beg) = lc(substr($$seq, $beg, $end - $beg));
  }

And @gaps is populated here:

	} elsif ($t[4] ne '.') { # an INDEL
	  push(@gaps, [$t[1], length($t[3])]);
	}

So around indels, the software will make a window of lowercase letters.

And if you have no coverage for a region, I'm not sure you necessarily want the software to just tell you that it's reference, because what if there's a deletion?

[ekg]

If you have a haploid sample, and only one sample, and the VCF reflects that fact (e.g. you called with freebayes --ploidy 1) then you can use vcfgeno2haplo (https://github.com/ekg/vcflib/blob/m...geno2haplo.cpp) with an arbitrary window size to get a consensus sequence. It will be recorded in your "ALT" in the only line of the output VCF. This works with indels.

If you have diploid samples, you'll have to phase them somehow. This method won't work if your samples aren't phased.

For coverage, you can get statistics in a number of ways. I'm not sure the best for your usage.