29
From the manual of Velvet:
it must be an odd number, to avoid palindromes. If you put in an even
number, Velvet will just decrement it and proceed.
the palindromes in biology are defined as reverse complementary sequences. The problem of palindromes is explained in this review:
Palindromes induce paths that fold back on themselves. At least one
...
14
It's difficult to get this to go massively quicker I think - as with this question working with large gzipped FASTQ files is mostly IO-bound. We could instead focus on making sure we are getting the right answer.
People deride them too often, but this is where a well-written parser is worth it's weight in gold. Heng Li gives us this FASTQ Parser in C.
I ...
13
To expand on the answer above, in case it isn't clear, we show:
Why palindromic sequences must be even in length
Why palindromic sequences induce self-loops in a de Bruijn graph
Why self loops in a de Bruijn graph are problematic
1. Palindromic sequence ⇒ sequence is of even length
Idea: in an odd-length k-mer, its middle nucleotide is 'flipped' in its ...
12
MultiQC can merge all your different reports into a single one.
Which could be useful once you manage to know which QC tools to use.
12
In any scenario where depth of coverage is an important factor, PCR duplicates erroneously inflate the coverage and, if not removed, can give the illusion of high confidence when it is not really there.
For example, consider the following hypothetical scenario.
*
...
11
tSNE often offers better visual representation (separation) on such complicated data than PCA. As Micheal pointed out, computing a tSNE embedding over 20.000 gene dimensions is computationally unfeasible, so a number of PCs are normally calculated and these are used as input for calculating the tSNE. They are used in tandem.
As for global vs. local, we are ...
10
I found a post useful for this topic.
It explains the difference of coverage and depth. It also has a useful explanation on how to calculate coverage and depth.
Here is a copy of what the link says just encase the post is removed:
Depth of coverage
How strong is a genome "covered" by sequenced fragments (short reads)?
Per-base coverage is ...
9
The Global Alliance for Genomics and Health has been working on the issue of representing sequencing data and metadata for storage and sharing for quite some time, though with mixed results. They do offer a model and API for storing NGS data in their GitHub repository, but it can be a bit of a pain to get a high-level view. I am not sure if any better ...
9
There are several questions in your post I'll try to answer each one:
Is there any way to calculate how deep the sequencing is ?
See gringer's answer. TLDR: The depth of the sequencing is how many times each position has been sequenced.
What should be the optimum depth to get reliable data ?
The optimum depth depends on what you want to do with that ...
9
Don't do it: Your FASTQ file is either malformed or your FASTQ record spans more than four lines, which is allowed in FASTQ. For a detailed description of what can go wrong in FASTQ parsing see for example http://biopython.org/DIST/docs/api/Bio.SeqIO.QualityIO-module.html#FastqGeneralIterator.
If the FASTQ is malformed, then you should really ask yourself ...
9
one-liner
Here's a gritty one-liner to count the number of reads in a region if you have just one region that you want to investigate. Change the 1 in ($4 >=1) and the 500 in ($4 <=500) to set your window. Change "hg19" to your target sequence.
Note, this one-liner does not double-count reads because of uniq.
samtools view file_sorted.bam | \
...
8
Sequencing depth is typically calculated as the number of total bases sequenced divided by the number of bases in the target genome. An Illumina sequencing run with 2x125 bp reads and 500 million read pairs sequenced would be a sequencing depth of about 40X (assuming my calculations are correct for a 3 billion base genome).
The sequencing depth depends on ...
8
If you don't mind a bit of manual counting, then samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam will produce output where you can count the bases for that position. You could, of course, use the command line to do most of that automatically:
samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam |
cut -f 5 | tr '[...
7
bwa mem is newer, faster, and [should be] more accurate, particularly for longer reads.
From the bwa man page (presumably in Heng Li's own words):
BWA is a software package for mapping low-divergent sequences against a large reference genome, such as the human genome. It consists of three algorithms: BWA-backtrack, BWA-SW and BWA-MEM. The first ...
7
We routinely run both FastQC and FastQ Screen on all of our raw sequencing reads. FastQ Screen is a tool for detecting cross-species contamination. MGA is another similar tool.
There are then lots of QC tools specific to different types of data, most of which run after alignment. For example RSeQC (RNA data), Qualimap and many many others. Without ...
7
Several papers have made this distinction, and a few indeed use different terms to distinguish between them. For example, Kazaux et al. (2016) acknowledge that:
These constraints favour the use of a version of the de Bruijn Graph (dBG) dedicated to genome assembly – a version which differs from the combinatorial structure invented by N.G. de Bruijn.
...
7
PCR polymerases introduce errors. When an error arises in the first few cycles of amplifications, it will appear in a reasonably high fraction of DNA fragments in the library. After sequencing, you may see the same error occur to multiple reads. If you remove PCR duplicates when calling variants, all errors are reduced down to one read. For high-coverage ...
6
I am working on a Illumina sequencing simulator for metagenomics: InSilicoSeq
It is still in alpha release and very experimental, but given a multi-fasta and an abundance file, it will generate reads from your input genomes with different coverages.
From the documentation:
iss generate --genomes genomes.fasta --abundance abundance_file.txt \
--...
6
I think the best method or combination of methods will depend on aspects of the data that might vary from one dataset to another. E.g. the type, size, and frequency of structural variants, the number SNVs, the quality of the reference, contaminants or other issues (e.g. read quality, sequencing errors) etc.
For that reason, I'd take two approaches:
Try a ...
6
The following is more than twice as fast; however, wc counts newline characters as well. We thus need to subtract the line count from the base count:
fix_base_count() {
local counts=($(cat))
echo "${counts[0]} $((${counts[1]} - ${counts[0]}))"
}
gunzip -c "$file" \
| awk 'NR % 4 == 2' \
| wc -cl \
| fix_base_count
However, the ...
6
pigz | awk | wc is the fastest method
First off for benchmarks with FASTQ it's best to use a specific real-world example with a known answer. I've chosen this file:
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01815/sequence_read/ERR047740_1.filt.fastq.gz
as my test file, the correct answers being:
Number of reads: 67051220
Number of bases in ...
6
As @AaronBerlin mentioned, you didn't remove reads that were completely trimmed. Next time use the --minimum-length option and set it to something reasonable, like 20. Alternatively, use "Trim Galore!", which is a wrapper around cutadapt that has more reasonable defaults.
6
If you are 100% sure the read only has 4 lines (they can have more), you can use this sed command:
sed -i.bak '/^@HWI-D00466:116:CC62WANXX:3:1102:7363:63646 1:N:0:GCACACG/,+3d'
The -i.bak makes sed modify the original file and create a backup copy with the same name and the extension .bak. The command just means "delete the line matching the pattern and ...
6
If you're lucky the information is stored in the BAM header in the lines starting with @PG which stores the command lines run on the BAM, use samtools view -H your.bam to display the header.
5
The quality control of ngs reads is heavily dependent on type of the project.
For genome assembly projects based on short reads, beside already covered checking quality of sequencing, you would like to look at the kmer spectra to find out, if your reads are going to make sense when they will be translated to De Brujin graph. You will get also a clues about ...
5
I personally don't think BQSR has a huge impact on variant calling, but you don't really need to guess. If you run GATK BQSR, it outputs a table and charts of exactly how much quality scores are adjusted. The adjustment will vary depending on the position in the read and genomic context (previous and following base). In my experience, the difference is a few ...
5
For metadata, I would use a SQL schema something like the following:
CREATE TABLE Project (
ac TEXT, -- project/Study accession
PRIMARY KEY (ac)
);
CREATE TABLE Sample ( -- biological sample/biopsy
ac TEXT,
PRIMARY KEY (ac)
);
CREATE TABLE AnalysisSample (
prj_ac TEXT, -- project acccession (Project.ac)
symbol TEXT, -- a short name ...
5
I agree that there is no ideal data model that is going to be stable for very long in a quick-moving field like genome informatics. Perhaps a schema-less (NoSQL or some other document-based system, such as MongoDB) database approach would work better? This gives you ultimate flexibility to attach whatever information is relevant to database entries you're ...
5
While your question is specific to cancerous germline mutations, I'd suggest you look at the COSMIC database of somatic mutations to include in your analysis.
There are other factors to include in this kind of analysis you're suggesting, such as predictive deleterious effects (PolyPhen for example can perform such predictions).
If you have 10M variants/...
5
Here are my attempts at definitions:
Sanger: A method of sequencing that depends on chain-terminatiing dideoxynucleotides. This sequencing uses the differential flow of DNA sequences of different lengths through a gel to determine the original DNA sequence, producing a single sequence per reaction container.
NGS: Next-generation sequencing, also referred ...
Only top voted, non community-wiki answers of a minimum length are eligible