Bcftools has sample/individual filtering as an option for most of the commands. You can subset individuals by using the -s or -S option:
-s, --samples [^]LIST
Comma-separated list of samples to include or exclude if prefixed with "^". Note that in general tags such as INFO/AC, INFO/AN, etc are not updated to correspond to the subset samples. bcftools ...
Most of these use RNA-seq data, some use WGS data, and some use both. They are listed alphabetically. I will add to the list when I discover more.
If you just want to filter out calls present in dbSNP then use:
java -jar GenomeAnalysisTK.jar \
-T SelectVariants \
-R reference.fasta \
-V patient.vcf \
--discordance dbSNP.vcf \
--discordance will produce calls not present in dbSNP.
Use VEP, ANNOVAR or snpEff to annotate your VCF file (I'd recommend combining your VCFs into a single file if they're single sample VCFs or are all comprised of samples from the same experiment/cohort). You can annotate the VCF file with 1000g (among a ton of other annotation sources). Once done, you can use bcftools view to subset the VCF as required.
They created 10kb bins, which could be called "dividing the genome". These are just 10kb long contiguous stretches (so position 1-10000,then 10001-20000, etc.). This is, of course, a nonsensical way of going about things, because "counting those up" in the regions you posted is the same as taking the width of the regions, dividing by 10000, and rounding a ...
If you have issue with memory and dealing with large object, maybe data.table is the way to go (https://github.com/Rdatatable/data.table/wiki):
Let's take two dummy data.frame:
dfA <- data.frame(Chrom = rep("chr1",150),
Pos = 1:150,
ref = rep("A",150))
dfB <- data.frame(Chrom = rep("chr1",10),
There is like a thousand different ways how to achieve this. You could use a specialised software for this (like bedtools) or calculate it simply in R.
R solution: You can make a function that calculates the number of SNPs in a range and than you can apply it on all the ranges in the table with genomic ranges.
snp_table <- read.table('WTSI-OESO_121_1pre....
You received an answer from Kamil for a similar question here: https://bioinformatics.stackexchange.com/a/8532/650
You can adapt the code they gave you to compare Chrom1 == sv & Start1 >= sv & End1<=sv to find events that happen within a gene. For events spanning genes, you're going to need to write a collapsing function for the gene ...
I have previously estimated tumour purity with the EXPANDS an inferred tumour heterogeneity program which is designed to calculate the number of clonal subpopulations in matched tumour/normal samples. The purity is essentially the size of the largest subpopulation identified in that sample - this is discussed in the programs FAQ. In addition to a matched ...
Your approach is not recommended. Always align against the entire genome. If you align against a subset you might accumulate false mappings. The aligner will always try to find the best match for every read. If the true origin of the read is not in the reference then the aligner will still try to find an acceptable mapping position. Therefore, align against ...
"DNA data" comes in several forms dependent on the technology used to produce them.
Companies like 23andMe are using SNP chips and those are available for an only rather a limited number of species. To be completely honest, I don't even know what would be the raw data they would provide, but I would suspect it could be a vcf file.
If the chosen ...
It's usually CNV callers that make use of Tumour/Normal WGS pairs to estimate purity. It can also be done with WES (exome) Tumour/Normal pairs.
There are several tools out there, I have some experience with the one written by Illumina (public on Github):
It requires realigning things with bowtie2, so I don't think it can ...
If you have gVCFs, the first thing you should try is joint variant calling. According to GATK, joint variant calling "empowers variant discovery by providing the ability to leverage population-wide information from a cohort of multiple sample[sic], allowing us to detect variants with great sensitivity and genotype samples as accurately as possible." source.
I'm not familiar with the program, but apparently Hail is setting itself up as a swiss-army chainsaw project for doing downstream analysis on variant-called datasets.
An overview of Hail can be found here:
A tutorial on association testing can be found here:
UPDATE: As an update, Sarah Walker (co-author on the poster) responded to my question on the GATK forum. She clarified with the following statement:
We believe the sites around 30X (and above 150X) are being filtered
due to low mapping quality (since it is whole genome so there are many
areas that are hard to map), which explains the low ...
In addition to the answer from @gringer there is a bcftools plugin called split that can do this, but gives you the added ability to output single-sample VCFs by specifying a filename for each sample.
$ bcftools +split
About: Split VCF by sample, creating single-sample VCFs.
Usage: bcftools +split [Options]
-e, --exclude EXPR ...
Thanks Christopher Chang for the great answer via the plink2 google group! (See here)
His answer was as follows:
"Yes, plink 1.9 did change the default merged sample order. However, you can request plink 1.07’s behavior by adding “--indiv-sort 0” during the merge. (You can also use --indiv-sort + --make-bed at any time to change the sample order to ...
You can't know which is the driver and which is the carrier. At most you can say that a specific gene deviate more of the expected underlying hypothesis. See also other resources online.
You also seem to ignore other genes that deviate more of your null hypothesis. I recommend to plot the histogram of the p-values to see if the distribution is uniform
I would suggest you look at the website for the excellent MultiQC tool (https://multiqc.info/). We use it in a clinical setting to collate all the QC statistics from our pipelines into a single report (It displays results from multiple QC programs such as FastQC, Picard, bcl2fastq, VerifyBAMid etc). The homepage has a selection of example reports - ...
If we save the script pasted in the main post as a sh file and we have some .vcf files in a folder, by this line we can iterate over vcf files to extract what mentioned in the script returning .txt output by the name of corresponding vcf
for file in *.vcf ; do bash indel_vcf_parasing.sh $file ; done
What is the difference of ReadCount in INFO column and DP in FORMAT column?
I'm not super familiar with Strelka but I believe it follows the VCF spec, so DP in the FORMAT column should give the depth of coverage (ie read count) per sample. ReadCount in INFO is not in the VCF spec but based on the definition it sounds like the depth of coverage aggregated ...
To help us help you, please make your data available with dput() next time.
To achieve your goal with ggplot2, you would need all of your data in one data frame and in the "long format".
dat1 <- "
Patient DEL INS SNP MNP total
LP6008337-DNA_H06 927 773 40756 0 42456
LP6008334-DNA_D02 1049 799 31009 0 32857
dat2 <- "
Patient DEL INS SNP MNP total
There seem to be WGS assemblies that can't be found in NCBI Assembly database, e.g.: https://www.ncbi.nlm.nih.gov/nuccore/1779902990. I guess that such assemblies also cannot be found in The assembly_summary.txt files that are described in https://ftp.ncbi.nlm.nih.gov/genomes/README_assembly_summary.txt.
My current best guess is that each WGS assembly has a &...
Are you looking too far ahead into the future? No. This is certainly possible now. The gold standard being Illumina or BGI short read WGS. Long read sequencing can capture some extra data but is very noisy.
How or where do one obtain such DNA data? There are now many companies. You want one the uses NGS with reads at least 150 bp and garantees at least a ...
You are correct that differences in ploidy between autosomes (always ploidy=2) and X chromosome mean that the same VQSR model cannot be applied to both cases equally.
However, this is the case only for male X chromosomes (in human species).
Therefore, yes, I'd recommend building a separate VQSR model when it comes to male X chromosomes. Importantly, the ...