Hot answers tagged

8

Could you please show us the context in which this appears, as you seem to be interpreting this differently to Devon. If it's appearing as you say, GA and AG, then Devon is right, this usually means that the sequence goes ###AG### or ###GA###, which are two very different sequences. If, however, as you're implying, the sequences are actually ###G### or ###...


8

I wonder if there is a simpler solution recently? (and hopefully, I can solve it within the scope of python. ) A simpler solution, I don't know... but this is at least one Python solution using Biopython's ELink method via NCBI's Entrez E-utils. The Biopython library is flexible enough (they have an in-depth tutorial worth reading) to modify the code below ...


7

I can show you a simple way in R, using biomaRt. Let's say you have two snp ids you want to exmine. SNPids <- c("rs431905509", "rs431905511") library("biomaRt") # To show which marts are available listMarts() # You need the SNP mart mart <- useMart("ENSEMBL_MART_SNP") # Find homo sapiens listDatasets(mart) # This will be the dataset we want to ...


6

What you are attempting to do is known as LD-pruning. As @Emily_Ensembl said, it is not customary to do this for standard association tests: it is possible that one of the SNPs you remove is causal, or a better proxy for the causal locus, and would give (slightly) better association signal than the other. Even for SNPs in perfect LD, pruning is unwise ...


6

Simulating genotypes with realistic correlation structures is indeed not so simple, and there's quite a few papers dedicated entirely to that (e.g. https://bmcgenet.biomedcentral.com/articles/10.1186/s12863-015-0173-4). Also, DEPICT (https://data.broadinstitute.org/mpg/depict/index.html) comes with a number of simulated GWASs to generate the nulls, so that's ...


6

You just searched for X[Chromosome], you didn't specify a species. Presumably, your dataset comes from a specific species, so you should limit your search to that species. For instance, for human, you can search dbSNP for X[Chromosome] human[organism]. Then, click on the send to link at the top left and download the data: This is veeeeeery slow, however. I ...


6

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 \ -o patient.filtered.vcf --discordance will produce calls not present in dbSNP.


5

What you are looking for is SNP annotation. If you have the chromosome:position reference and alternate alleles for your SNPs of interest, it can be as simple as uploading them to the variant effect predictor. http://grch37.ensembl.org/Homo_sapiens/Tools/VEP This will give you the predicted protein change and novelty of the variant with respect to known ...


5

Most likely, the GWASs that generated your summary statistics used other imputation panels than 1000G, like HRC. Clearly, PLINK can't estimate the LD for SNPs that aren't found in the reference, and most likely just ignores them. I think the easiest solution is to use whatever individual-level dataset you have to estimate the correlations - even if it's ...


5

I know that the GWAS association p-value threshold is 1e-8 This may be a common threshold of statistical significance that is used, but it's definitely not an absolute value. It's a hack to try to work around many issues with GWAS associated with testing millions of SNPs. Unfortunately, the most relevant issue in GWAS (for spurious significance) is ...


5

(edit)you can filter the VCF annotations with snpsift, I've also written a VcfFilterSequenceOntology http://lindenb.github.io/jvarkit/VcfFilterSequenceOntology.html I've written vcf2table: http://lindenb.github.io/jvarkit/VcfToTable.html It decodes VEP and SNPeff annotations: >>chr1/10001/T (n 1) Variant +--------+--------------------+ | Key | ...


5

In theory, almost any base in the human genome may mutate, so you have billions of variants to go. Ok, this is not so useful. A related and potentially useful question is: given a human, what is the fraction of his/her variants seen in the the previously sequenced samples? For a Wright-Fisher population, there is an analytical answer: x% of variants on a ...


5

As far as I'm aware, Illumina provide CSV annotation files for all their sequencing chips, which can be used when they can't be found in Bioconductor. You can find annotation information for the PorcineSNP60 here, in particular the Manifest file (CSV format). The format is Illumina's weird "we say it's a CSV because there are commas in it" format, so if ...


5

I’m no longer working in tumour sequencing so I’m by no means an expert. But in a nutshell, the reason is that, as indicated, homozygous SNPs aren’t informative: if your allelic fraction is 100%, meaning you only observe a single nucleotide at a given position, we don’t know whether we’re dealing with a tumour related mutation. In fact, we probably don’t ...


5

The difference between the two depends on to whom you talk ;) You are right: both refer to one base difference from the sequence. Usually the term "mutation" is used if the change has an impact on the phenotype. The "P" in "SNP" means "polymorphism". A lot of people use this term if more than 1% of the individual in a given population have this variant. My ...


5

Just use bcftools view for filtering: $ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz > output.vcf To truncate this list to 65,000 SNPs count the header lines, sum them to the number of SNPs you like to have and use head -n. $ bcftools view -h input.vcf.gz|wc -l 255 $ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz | ...


4

Try this pipeline: Clean raw fastq files (trimmomatic) Align clean fastq files to the reference genome (dna-to-dna aligner of your choice) Convert the alignment sam file to bam (samtools) Call SNPs from bam file (samtools) Here's a start tutorial on this: http://ged.msu.edu/angus/tutorials-2013/snp_tutorial.html


4

Even with inbreeding and other genetic phenomena that might mask actual evolution of these cultivars, any phylogenetic methodology would be capable of determining relationships accurately. Try creating a Neighbour Joining tree with MEGA, which is one of the simplest methods available. This should give you enough to check the relationships of the cultivars.


4

Is "user5054" the same as "user5504"? No? Exactly. Not only does order matter, it's incredibly vitally important. AG and GA are completely and totally different from each other. If this is in a coding region, then the resulting amino acid is undoubtedly changed (fun fact, the only exception is TAG and TGA). If it's at a splice site then splicing is likely ...


4

You can assume that the overwhelming majority of rsIDs are the same between GRCh37 and GRCh38 (they're semi-stable IDs). There are, however, a number of rsIDs that are present only in GRCh37, which you can find here. Note that the format of this file is a bit strange, it's chromosome|position|ID|weight, where position is sometimes empty, weight is typically ...


4

If you mean that you want snps from individuals, instead of all together, you can find them in 1000 genomes. Here different individuals from different populations are sequenced and variants are called, including snps. Maybe you should add to your question what your goal is, what would you like to find?


4

OK, I ended up contacting one of the authors of the program, who kindly answered to my question: Either approach should work, [i.e. running the analysis by chromosome vs joining the data together] but I suggest dividing by chromosome because that approach has had more testing. I hope this might help other users, I saw that several questions ...


4

snippy-core is a tool that processes the snippy_outdir_{1..n} folders and produces the following files: core.aln – a fasta format file that contains only the polymorphic sites (of {'A', 'C', 'T' and 'G'}) in your input set, whose length tends to decrease asymptotically as more isolates are added core.full.aln – a fasta format file that contains all the core ...


4

From my memory of what a statistician told me, a PCA aims to determine independent linear combinations of variables (i.e. genotypes) that account for the most variation in the dataset. With 10 million SNPs, the vectors describing the genotypes of all individuals reside somewhere inside a 10-million-dimension hypersphere, and the first principle component ...


4

The simplest approach is to remove the first occurrence of a _ followed by a capital letter, another _ and another capital letter from each line: $ sed 's/_[A-Z]_[A-Z]//' file 1 rs2880024 866893 A G 1 chr1_867635 867635 A G 1 chr1_869303 869303 A G 1 chr1_873558 873558 C A 1 ...


3

ANNOVAR is another tool that will help you functionally annotate genetic variants. It will tell you if the variant is known in COSMIC or Clinvar databases, which AA change will occur, if it's a frameshift, synonymous, non-synonymous, etc...


3

You can use imputation for guessing at what the missing SNPs might be based on known LD patterns in populations. This procedure will give you an idea of whether the recombinational history of genomic regions is important, rather than specific SNPs. Given that you're looking at LD-based measures, this should be okay (although the logic behind the calculation ...


3

If you only have one gene and you only need to do this once then the simplest possible workflow is to generate the aliment using STAR (optimally with the two pass method) and open the two resultant .bam files in IGV, the coverage column at the stop above your alignment track should clearly show the counts of reference and non-reference supporting reads. The ...


3

For counting reads I use mpileup, e.g. samtools mpileup --reference hg38.fa -r Chr10:18000-45500 input.bam, which will give base-resolution coverage for a BAM file. I've written my own script to process mpileup output and make it easier to understand. By default it reports read coverage as a proportion of total coverage, but this can be modified by using ...


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 ...


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