23

You're the second person I have ever seen using NCBI "chromosome names" (they're more like supercontig IDs). Normally I would point you to a resource providing mappings between chromosome names, but since no one has added NCBI names (yet, maybe I'll add them now) you're currently out of luck there. Anyway, the quickest way to do what you want is to samtools ...


11

It’s a matter of preference I guess but I recommend the Ensembl builds. Decide whether you want the toplevel or primary assembly, and whether you want soft-masked, repeat-masked or unmasked files. The naming schema is very straightforward; the combinations are described in the README file, and all files reside in one directory. For example, if you want the ...


9

tl;dr: Just use the either the downloads on the Bowtie2 homepage or the Illumina iGenomes. Or just uncompress and concatenate the FASTA files found on UCSC goldenpath and then build the index. A bit longer answer: There are two components to "genome for a read mapper" such as Bowtie or BWA. First, you need to choose the actual sequence (genome release ...


7

GATK has a solution that might work for you: FastaAlternateReferenceMaker, which : "Given a variant callset, this tool replaces the reference bases at variation sites with the bases supplied in the corresponding callset records." Input The reference, requested intervals, and any number of variant ROD files. Output A FASTA file representing ...


7

tl;dr: technical difficulties, or sex For the gorilla, that genome is from a female gorilla (Kamilah the gorilla) so she doesn't have a Y chromosome. For the chimp link, there is a visualization of the chromosomes that you can click on that takes you to the genome browser for the Y chromosome (see image). Notably, there is a more recent genome build that ...


6

You'll get the exact same index (the amb, ann, bwt, pac and sa files) whether the reference is gzipped or not. BWA also makes its own packed reference sequence (the .pac file) so you don't even need the genome around after you index.


6

At the moment, the standard reference genomes (e.g. hg19, hg38) are haploid genomes. We know that the human genome is diploid. Naturally, the latter would be the respectively correct representation of the human genome. The premise of the OP's question is false. The natural reference representation of the human genome is not diploid. Think of a reference ...


5

For calling small variants, the standard way is to simply call diploid genotypes. You can already do a variety of research with unphased genotypes. You may further phase genotypes with imputation, pedigree or with long reads/linked reads, but not many are doing this because phasing is more difficult, may add cost and may not always give you new insight into ...


5

I will answer one of the points – E.coli. TL;DR Bacteria, and in particular E.coli, are highly variable and there is usually no single best assembly. Large scale WGS studies should come with multiple assemblies for individual monophyletic clusters. Long answer: Whereas a single reference sequence can make sense for human or mice (to a certain extent), ...


4

For a quick (but reliable) analysis, I'd recommend using Kallisto or Salmon to quantify isoform read counts using the transcriptome of the lab strain. If you have a concern about transcripts that are in your sample but not in the lab strain, you can do two Trinity assemblies: one fully de-novo, and one genome-guided. These transcripts can then be used as a ...


4

The "right" solution would be realignment, but that's expensive and most of us would not go that route. My preferred solution would be to convert the bed file, as opposed to the bam. Here's why: 1) Reheadering the bam means that you may have reads aligned to contigs without a corresponding entry in UCSC (see Devon's list for the mappings). This is a problem ...


4

You could convert VCF to BED via vcf2bed --snvs, vcf2bed --insertions, and vcf2bed --deletions, and then use samtools faidx by way of a wrapper script to convert BED to FASTA, e.g.: $ vcf2bed --snvs < variants.vcf | bed2faidxsta.pl > snvs.fa $ vcf2bed --insertions < variants.vcf | bed2faidxsta.pl > insertions.fa $ vcf2bed --deletions < ...


4

Your problem is caused by using the transcriptome fasta file rather than the genome fasta file. You've already given it transcriptome information with genome_genes.filtered.gtf, it needs the genomic sequence to go along with that. As an aside, your GTF filtering command is unlikely to be useful unless you have a Gencode or Ensembl GTF for human. Have a look ...


4

The output of gffcompare includes several files per run (just like cuffcompare). Example for a run: $ ls cuffcmp.* | sed 's/\t/\n/' cuffcmp.combined.gtf cuffcmp.loci cuffcmp.output.gtf.refmap cuffcmp.output.gtf.tmap cuffcmp.stats cuffcmp.tracking From my experience, the easiest file to manipulate is the tmap file: Tab delimited file lists the most ...


4

Yes, there are. There are some suggestions in the comments (VG, WhatsHap, GraphAligner, Minigraph). However, to be clear, the current default variant calling algorithm used in the GATK (the haplotype caller) is a graph-based algorithm. It constructs a De Bruijn graph from mapped reads for local assembly around variants. Other tools that also follow this ...


3

There's a vcf2fq sub-program that was written as part of vcfutils to convert a VCF file into a fastq file given a reference sequence. Unfortunately this doesn't work properly with INDELs (it will just mask them, rather than actually converting them), so I wrote a modification to implement INDEL correction as well: ./vcf2fq.pl -f <input.fasta> <all-...


3

As others have said, doing whole-genome (or even whole-chromosome) alignments is the wrong solution. Simply track where you're creating SNVs. If you wrote your SNV locations to a sorted BED file you could use something like the following: #!/usr/bin/env python minWidth = 100000 # Only report regions of >=100kb last = [None, None] for line in open("some ...


3

This seems relatively complicated given the structure of a BSGenome object. The creator of the package answered this question previously on the Bioconductor support forums: https://support.bioconductor.org/p/86665/#86757 We don't provide an easy way to inject arbitrary SNPs in an arbitrary BSgenome at the moment. However, it should not be too hard to ...


3

I couldn't find how is this package build, my guess is that is as in the example of the OrganismDBI package. See that vignette section to build your own data package with the genome you want. The locations are different because you ask for the transcript start and end, if you add the "TXNAME", you will (probably) see there are different transcripts, so ...


3

There are multiple ways. Here are two: samtools idxstats ‹bamfile› samtools view -H ‹bamfile› | grep '^@SQ' Both commands give you a list of the references in the BAM file. To get their number, simply append | wc -l. But note that idxstats also outputs a final entry, *, which doesn’t correspond to any reference but to all unmapped reads. So subtract that ...


3

First, to answer your question about mapping to a low-quality reference: 1. Mapping For mapping, low genome contiguity (low N50) doesn't really matter. You will be using a spliced aligner and short reads, so even small portions of your reads will map if they match the reference genome. What does matter though is that you are confident that the draft genome ...


3

I haven't done any whole-genome STR analysis from NGS data myself, but are aware of others that have used lobSTR for this. There's also a recent paper [here] that compares a few different STR analysis packages (i.e. RepeatSeq, LobSTR, HipSTR, GangSTR). Here's the concluding paragraph: In conclusion, all these tools are built to genotype STRs but have ...


2

The current ensembl entry doesn't have a 29 either. The archived ensembl assembly lacks 29 30, and 31 and 33 and LGE64. The chromosomes after 30 are tiny, so they might not be visible in a karyotype. They probably realized that "chr 29" was really attached to some other chromosome.


2

If you store the positions of your random SNPs in a sort-bed sorted BED file, you can filter with BEDOPS bedmap: $ bedmap --count --echo --range 50000 --delim '\t' SNPs.bed | awk '$1==1' | cut -f2- > SNPs.filtered.bed The file SNPs.filtered.bed will only contain those elements from SNPs.bed which do not overlap within 50kb. Or you can use BEDOPS bedops ...


2

First of all, take these data with a very large pinch of salt. This sort of targeted analysis is not designed to produce high quality genetic data but, to give an idea of a sample's ancestry. Given a file like the one you describe with no quality information, you have no way of knowing how reliable any of the variants called actually is. From personal ...


2

It's pretty common for SNPs and indels that map to the same position in the genome to be assigned different rsIDs, particularly when the indel spans multiple base pairs in the reference. You can get more info about these particular variants here and here. Or you can see how both of these variants relate to the reference sequence here. Without info about the ...


2

It sounds like you want samtools faidx foo.fa followed by samtools faidx foo.fa chr1 > your_subset_file.fa (or whatever the first chromosome is). The output file is then a regular fasta file subset as you like. You can get any chromosome you want with that. In fact, you can also do regions (e.g., chr1:100-1000, though note that the sequence name in the ...


2

Disclaimer, I am not an expert on admixture, answer is provided without any warranty. It depends what you want to figure out. Why do you want to categorize the ancestry of the individual? If you analyze this origin of population X you should make a reference without population X and analyze all the individuals from the population with all the surrounding ...


2

Like gringer already commented, the cause of your error is nicely given by the java program. Your default max heap settings are probably too low for the program. If you have more memory on your machine you can give max heap a higher value. See this post and answers on Stack overflow. You'll probably need to use the -Xmxn argument in java.


2

You should be able to use something like Pilon to convert from one genome to another, assuming the changes are all fairly small local changes (but that's almost never the case). I find your resistance to long-read data interesting. It's important to consider the cost of bioinformatics as well. Sequencing a phage genome with long reads will take less than an ...


Only top voted, non community-wiki answers of a minimum length are eligible