3
$\begingroup$

I'm having problems properly looking for a solution because I'm a layman in Bioinformatics not familiar with the terminology. I'm hoping you can nudge me in the right direction, please! Thank you very much.

What I want in the end

A FASTA file containing a single sequential sequence (chromosome 1-22..) of the human genome.

What I have

  • 23gb of "clean_data" which are .fq.gz files.
  • 21gb of "result_alignment" which is .bam+.bai
  • "result_variation" indel.vcf.gz and snp.vcf.gz

What I tried

I used samtools fasta reads.bam > reads.fasta which results in many lines like so:

>V100005418L4C001R001322578/1
AAATAGCCCACAC...
>V100005418L4C001R006084507/1
GCCATAAAGCCTA...
...

But these include many many duplicate lines and the reference genome from the human genome project GRCh38 is only 3GB in size, so I thought that this FASTA sequence needs to be properly "aligned" somehow?

Am I thinking correctly here?

I'm really confused with all of the terminology and don't know exactly what I need to look for. Could someone kindly put the process into simpler words for me?

Improve this question
$\endgroup$
3
  • $\begingroup$ What you have in BAM format is an alignment of reads to a reference. What you are looking for (a single fasta per chromosome) is a new assembly. Using "samtools fasta" will just get you each read in fasta format, which is clearly not what you want. In addition to doing a (de novo) assembly of your reads you could make a (reference-guided) consensus sequence. You can find some pointers on how to do that here: biostars.org/p/367960 $\endgroup$
    – Wouter De Coster
    Commented Jul 5, 2019 at 5:21
  • $\begingroup$ Thank you so much for your comment @Wouter. The link you posted sounds very promising. I will try it out and get back to you. Meanwhile I clarified in my post what kind of data I have. Maybe you want to take a deeper look. Thanks again! $\endgroup$
    – Daniel Biegler
    Commented Jul 5, 2019 at 14:05
  • $\begingroup$ another option is to take the reference genome you used for alignment, and 'mutate' it with the variants you identified. One tool to do that is GATK FastaAlternateReferenceMaker. I don't know how well it deals with indels though. And keep in mind that the sample which you sequenced is diploid. If you want one fasta sequence per chromosome for that individual then you have to think about how to handle heterozygous sites. If you want one paternal and one maternal fasta then you are in bigger trouble because you need to know on which phase the variants are relative to each other. $\endgroup$
    – Wouter De Coster
    Commented Jul 5, 2019 at 14:23

2 Answers 2

Reset to default
3
$\begingroup$

If you just want the human chromosomes in a FASTA format, why not downloading directly the individual chromosome files (chr*.fa.gz)?

If you don't rely on the "official" releases (for whatever reason), then what you need is first to assemble your reads into genomes, possibly (but not surely) going down to a chromosome-level resolution. Different approaches for assembly are possible/available, mainly depending on (1) what do you biologically need to answer and (2) the type of sequencing data you have - we'll need more details to be able to help you in this direction.

The samtools fasta utility will "just" convert your reads (likely given from a sequencing facility) from an alignment (BAM) format into a reads (FASTA) format, there is no concept of chromosomes here. You won't get assembled chromosomes this way, just extract the sequence in FASTA format from your alignment file.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you very much for your input! I don't have something biologically to answer, this project revolves about sequencing(?) this specific individuals DNA. I need the DNA from "start to finish" so to speak, Chr1 to Chr22, X,Y in ascending order, is what I'm trying to say. I understand now that the FASTA conversion is not what I need, I need to assemble my reads first. As for data I have: 23gb of "clean_data" which are .fq.gz files. 21gb of "result_alignment" which is .bam+.bai and I got "result_variation" indel.vcf.gz and snp.vcf.gz. Can I use these to achieve what I want? Thank you very much! $\endgroup$
    – Daniel Biegler
    Commented Jul 5, 2019 at 13:57
  • $\begingroup$ Usually when you sequence a person's genome, you're trying to identify variants that the individual has, compared to the reference genome. In that case, you should do what @aechchiki said, and download the official reference FASTAs, and create a VCF file containing all the variants for that individual. Alternatively, it also sounds like you might be trying to do de novo genome assembly. That is a very different topic from variant calling, and that would give you a "reference" that is very different from the reference genome builds mentioned above $\endgroup$
    – James Hawley
    Commented Jul 5, 2019 at 15:43
  • 1
    $\begingroup$ @DanielBiegler please edit your question and tell us why you are sequencing this individual since the data you will require will depend on what you are trying to do. It is extremely unlikely that you need the individual's full genome rather than the specific locations where their genome differs from the reference (variants). But if you don't us what biological question you are trying to answer using the individual's sequence, then we can't really help much. $\endgroup$
    – terdon ♦
    Commented Jul 5, 2019 at 17:10
2
$\begingroup$

What I wanted to do is called "Consensus Sequence". Two steps were needed:

$ bcftools mpileup -Ou -f ref.fa input.bam | bcftools call -Ou -mv | bcftools norm -f ref.fa Oz -o output.vcf.gz

After that you can create the consensus:

$ bcftools consensus -f ref.fa output.vcf.gz > out.fa
Improve this answer
$\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.