cannin/biostars_qa huggingface dataset:
Question: Hi.

I am wondering a term with SIFT score.

I think that SIFT refers to some measurement of SNPs, and while reading Annovar paper, I saw below sentence as follows:

> Finally, Annovar can filter specific variants such as SNPs with >1% frequency in the 1000 Genomes Projects, or non-synonymous SNPs with SIFT scores > 0.05.

Regarding above sentence, I ask you two questions.

1. I think that 1% frequency is a little bit low allele frequency. Dose it have an effect to filtering irrelevant snp variants? I don't think so..
2. SIFT-score threshold is about 0.05 as shown in above sentence. What does SIFT means about and threshold of 0.05 might be effect on filtering variants?

Answer: SIFT and PolyPhen are the two most commonly used algorithms for predicting if a SNP has a (generally negative) effect on protein structure. Due to the nature of the redundant genomic code, many SNPs never translate into any effect in the protein - far more than you would expect by chance - because variations which effect protein sequence are usually under negative selection pressure - so SIFT/PloyPhen can be used to weed out a lot of irrelevant stuff from a very large list of candidate variations.

If i'm not mistaken (and its been a long time since i used either, so i might be making this up) SIFT's algorithm gives more weighting to variations which change the net charge of the protein, while PolyPhen uses aminoacid or base conservation to determine relative importance. Both obviously rank premature stop variants and other nonsense variants very highly - so often there is a lot of overlap.

Again, it's been a long time since I used either, i might have gotten that the wrong way around. But i can tell you this - I spent 3 years studying single-basepair exon variants in consanguineous families with a known phenotype, and very very very rarely did SIFT or PolyPhen ever guess the correctly the variant from a list. I wouldn't say they are junk, they're not, but variants which caused transcription factor non-specificity, splicing variants, RNAPol destabilisation, etc, are completely ignored. Do not rely on SIFT and Polyphen for anything other than ordering a candidate list for follow-up analysis :)
Question: I have multifasta sequence and i am trying to rename the header and print the first id only

>gene_1|id=abc|protein=capA|ABD5672.4
ATGCGTCGACGTCGTACGGGTTTT
CGTACGGGTTATGCGTCGACGTCG
TACGGGTTTT

desired output:

>gene_1
ATGCGTCGACGTCGTACGGGTTTT
CGTACGGGTTATGCGTCGACGTC
GTACGGGTTTT

I tried following using sed but didnt work it seems sed -r 's/^([^|]+\|){0}//;s/\|//'any improvement on this? or any other alternative ?

Answer: Is this what you want?

`cat gene.fasta | sed 's/|.*//'`

It'll just remove anything after the pipe symbol `|` which separates the field you want from everything else
Question: I am conducting a differential gene expression analysis to determine the difference in the change in gene expression in each sex due to a treatment in human samples. For this type of study would you recommend UCSC or ensembl? I am concerned by the significant difference in GO terms and kegg pathway results that I get from alignment and differential expression analysis using GRCh37 ensembl genome vs hg19 genome for human samples. There was a significant difference between the two all the way through the workflow, which was to be expected since the two reference genomes and annotations are different. I was wondering which one is most reliable? My results are so different between the two that I'm finding it hard to make sense of it all. Any help would be greatly appreciated.

Answer: Could you please be more specific with your workflow, tools used and input files. Maybe you can share some of the input files and intermediate results with us? Overal, huge difference in GO terms for same genes in different annotations is strange.

You obviously looked at [https://www.biostars.org/p/113100/][1] , right?

[1]: https://www.biostars.org/p/113100/
Question: I want to store a docker image with the input data, code, software, intermediate results and results. This will be a massive file.

Is there something like the sequence archive for people who want to upload a docker image of their whole study? If not, what is the best/cheapest way to store such a file?

Any related thoughts not directly answering the question are also welcome.

Answer: Storage of very large Docker images is potentially a bit of an issue. As Devon suggested you can check into Quay.io, but I don't know if they have any file size restrictions. In practice downloading a large single file can be a bit of a problem. You might consider separating the data storage from the docker image itself that replicates the code aspects and software. Then you could put the Docker image anywhere (Quay.io or just Docker Hub) and store all of the data in a public Amazon S3 bucket (or another cloud storage system with a nice API). The Docker Image could have code for pulling the data from AWS
Question: I have a fastq file and need to remove the identifier, quality, and white space. So far my script looks something like this:

#!/usr/bin/perl -w

use strict;

my $inputFileName = shift;
my $outputFileName = shift;

open INPUTFILE, "<$inputFileName" or die "poop";
open OUTPUT, ">$outputFileName" or die "poop";

my @bases = ('A', 'G', 'T', 'C');
my $line;

while ($line = ) {
chomp $line;
if ($line =~ /^\s*$/)
elsif ($line =~ /^\s*@/)
elsif ($line =~ /^+/)
else {print OUTPUT $line, "\n";
}
}

However I keep getting an empty output file. I'm very new to perl, so be gently.

Thanks!!

Answer: #!/usr/bin/env perl
use strict;
use warnings;

my $lines = 4;
my $delimiter = "\t";
my $input = shift @ARGV;

open(INPUT,"$input");
while () {
if($. % $lines == 2){ # the % character is the modulo operator Pierre was talking about
print
}
}
close(INPUT);

Put that in a file called `extract_reads_from_fastq.pl` and do the trick with:

perl extract_reads_from_fastq.pl input.fastq
Question: Hi

I have a folder with 4000 files in it. I would like make a new folder, and copy into it all of the files which contain >50 fasta sequences. How do I do this?

I know that I need to create a simple loop, then grep '>' | wc -l and select only those with >50. But I am new to programming and am unsure how to write this properly.

Answer: loop over the files in DIR1,
use awk to print the files having more than 50 fasta sequences,
read and copy those files

$ ls DIR1/*.fa | while read F; do awk '/^>/ {N++;} END {if(N>50) printf("%s\n",FILENAME);}' "${F}" | while read F2; do cp "$F2" DIR2 ; done ; done
Question: I made individual list file like this:

N2002620115
N2002774440
N2002570297
N2002487409
N2002348348

And used commands

--bfile data --keep ind_list.txt --make-bed --out

According to plink web site tutorial, it must have paired columns - Family ID and Individual ID.
I made one column file and two columns file but all failed.

Reading individuals to keep [ data/ind_list.txt ] ...
ERROR: Problem with line:
N2002620115

on my plink screen.

How can I fix it?

Answer: *ind_list.txt* should have 2 columns FID and IID (FamilyID, IndividualID).

https://www.cog-genomics.org/plink/1.9/filter

> **--keep** accepts a space/tab-delimited text file with family IDs in the first column and within-family IDs in the second column, and removes all unlisted samples from the current analysis. **--remove** does the same for all listed samples.

If we don't have FamilyIDs then just paste the file to make both IDs same:

paste ind_list.txt ind_list.txt > ind_list_new.txt
Question: Hi,

How should I download the [dbGaP CF Genotype-Tissue Expression Project (GTEx)][1] dataset, it looks a very interesting data sets contains a variety of phenotypes and tissue attributes.

I am looking forward to do research on it, but I could not figure how to download it in the first place. Could someone please help me?

[1]: http://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/dataset.cgi?study_id=phs000424.v6.p1&phv=169063&phd=3910&pha=&pht=2743&phvf=&phdf=&phaf=&phtf=1&phs_code=&dssp=1&consent=&temp=1

Answer: From your link, click the "Study" tab, Jump to "Authorized Access" and follow the "Data access provided by:" link. You will be redirected to the following page, that has all instructions on how and who can apply:

https://dbgap.ncbi.nlm.nih.gov/aa/wga.cgi?page=login
Question: Hi, apologies for this basic question as I am new to the field. I have been checking my NGS data using FastQC and the checks fail on the Kmer content section. There appears to be a sequence TAGATCGGAA at position 90-100 bp in the reads which is enriched around 12 fold (Obs/Exp Max). However nothing shows up in the 'Overrepresented sequences' or 'Adapter content' sections of the report (complete flat line). The sequencing was done as BGI and I do not know what primers were used. My question is should I be trying to remove this Kmer sequence?

If I use grep on the first million reads in my fastq file to look for the sequence I only find 36 which appears to be quite a low number?

$ gunzip -c data1.fq.gz | head -4000000 | grep TAGATCGGAA | wc -l

Thanks for the help

Answer: No, there's no point in removing that. It's likely that that's just a bit of adapter contamination (I think the adapter contamination section of FastQC is looking for closer to full-length adapters, rather than just 5-8 bases). Go ahead and trim adapters regardless.
Question: I invested some time into learning R and would prefer not to learn Python or Ruby just for convenience of parsing XML.

I am not a pro R programmer so I need a higher level library then XML package of R.
getNodeSet and XPath is very hard to use.

Answer: One option is rentrez, from the rOpenSci family of packages.
Question: Hi,

I am trying to build a pipeline for detecting de novo mutations for trios WES data.
After obtaining the VCF file, and performing all the possible recalibration and refinement steps, I reached this step where I have to phase my variants.

As I have understood, please correct me if the following is wrong:

1- phasing means relating the genotypes to their paternal and maternal origin

2- The de novo mutations will not be phased, because they violate the mendelian law of inheritance

3- PhaseByTransmission can be used to do this job perfectly: i.e. to phase varinats and find denovos

Thanks in advance

Answer: 1) There are two common types of phasing: phasing by pedigree (by relating the alleles to the paternal or maternal origin via allele transmission), or local phasing (by determining which variants are locally in phase with each other, either by something like read backed phasing, or variant callers that directly call local haplotypes).

2) The de novo mutations can not be directly phased by pedigree, but they can be locally phased with respect to nearby variants (which may themselves be phased by pedigree).

3) You could use PhaseByTransmission, but you will obtain a better overall result if you jointly call the family using a pedigree-aware variant caller, such as `rtg family` or (for larger pedigrees) `rtg population` from [RTG Core][1]. This is because the pedigree-aware joint calling allows the evidence for each of the samples to influence the calls in other members of the pedigree (in `rtg population`, you can even use this to impute genotypes for missing family members during calling, which gets better the more family members you have). These callers automatically phase the output variants according to the pedigree, and directly output VCF annotations indicating putative de novo variants (including a de novo specific score).

(RTG Core also includes the `rtg mendelian` command which is useful for annotating VCFs for mendelian inheritance errors etc, and this command is also present in the smaller utility package [RTG Tools][2]). Disclaimer: I work for RTG :-).

[1]: https://github.com/RealTimeGenomics/rtg-core
[2]: https://github.com/RealTimeGenomics/rtg-tools
Question: Hi,

Would you please recommend me a good software for searching domains in a protein or coding nucleotide sequence? (except NCBI conserved domain search). Your help is greatly appreciated.

Answer: I decided to use Prosite for a project I'm working on because it has the best integration with UniProt. One of the problems with Pfam is that it is not updated nearly as often as UniProt is, and so your domain intervals are sometimes out of synch. Prosite is updated with UniProt, and you can pull the Prosite domains directly from UniProt, i.e.

http://www.uniprot.org/uniprot/?sort=score&desc=&query=id:Q12851&format=tab&columns=id,feature(DOMAIN%20EXTENT)

While I'm not aware of any out of the box packages for searching domains, you can write a custom script to pull the domain information you need fairly easily.
Question: Hello!

I've been trying to find a software that does what the post title describes, but without luck. I tried installing JalView but it won't install for some reason "Installer User Interface Mode Not Supported".

I have hundreds of images to generate so going to a browser MSA viewer and manually saving is not feasible. Am I out of luck or just bad at googling?

Big thanks in advance!

Answer: There’s a commandline version of `Espript` IIRC. That would probably work.

I also like the TeX package `texshade` for rendering MSAs, if you’re prepared to get your hands dirty with TeX.

I’ll post an example when I’m at a proper computer.

Finally at a computer, so here's an example of something you can do with TeXShade. You can of course run TeX from the commandline. It would probably be sufficient to write a wrapper script to inject the file names, or to actually write the source itself for you.

\documentclass[a4paper, oneside, 11pt]{report} % Doc Setup
\usepackage[left=1cm, right=1cm, top=1cm, bottom=1cm]{geometry} % Doc Geometry
\usepackage{texshade}
\usepackage{lscape}

\begin{document}
\begin{landscape}

\begin{texshade}{/path/to/your/alignment.ext}
\seqtype{P}
\setfont{numbering}{tt}{md}{sc}{small}
\setfont{names}{tt}{md}{up}{tiny}
\setfont{residues}{tt}{md}{sc}{tiny}
\setfont{ruler}{tt}{bf}{sc}{tiny}
\shadingmode[allmatchspecial]{similar}
\feature{top}{1}{25..25}{fill:$\downarrow$}{Hi, I'm an annotation}
\feature{bottom}{1}{70..76}{brace}{So am I!}
\end{texshade}
\end{landscape}

\end{document}

This creates:

It might not be a perfect solution, as you'd need to play with the sizing etc, or perhaps figure it out dynamically. I *think* `Espript` actually runs LaTeX behind the scenes anyway.
Question: Dear all,

during a step of a CNV analysis on cancer genomes, I would like to be able to COUNT and DISPLAY along the CHROMOSOME AXIS the NUMBER of READS in 10KB REGIONS of GERMLINE and CANCER GENOMES.

I would like to ask you please for your suggestions about :

-- any fast function/algorithm that COUNTS the reads from a BAM file in 10kb windows ?

-- a package/function that DISPLAY the COUNTS along the chromosome axes ?

many thanks,

bogdan

Answer: Have a look at [indexcov][1], but its intervals are of 16kb.

[1]: https://github.com/brentp/goleft/tree/master/indexcov
Question: Hello,

I have a bfile data that contains zero "0" in the allele column in .bim file.
Both allele columns contain SNPs recorded as "0". I didn't find corresponding descriptions on plink website. Does it mean unknown allele? Why is there such thing?

example:
> 1 rs123 0 242345 0 A

> 1 rs234 0 24532 0 0

Thanks!

Answer: In a .ped + .map fileset, it's possible for a variant to have only one distinct allele (everyone's A A or 0 0), or even zero (everyone's 0 0). If a .bim file is generated from this type of data, it'll have 0 in place of each unknown allele code.
Question: I have a very long list of coordinates that looks like that: `genome_build:chrN:coordinate`

I want to know if the given coordinate hits a gene or whatever else. What is the most efficient way to do so using python?

Example: having hg19:chr8:119331103 as input, it should return [SAMD12]

Answer: Finally, I wanna share the best solution I have found so far. 'Best' because it is extremely simple and pythonic.

All you need is [PyEnsembl][1]. Here is how easy it goes:

```
from pyensembl import EnsemblRelease

data = EnsemblRelease(75, auto_download=True) # will take a while to get the genome data installed in your system
gene_names = data.gene_names_at_locus(contig=6, position=29945884)
print gene_names

>>[u'AC023590.1', u'SAMD12']
```

[1]: https://pypi.python.org/pypi/pyensembl/0.5.4
Question: Hello everyone,

I've looked at `vcftools`, `bcftools`, `bedops`, `bedtools` and GATK `SelectVariants`, but I don't see how I can address my requirement - I hope the community can help.

I want to subset a VCF file based on 4 attributes: `CHROM`, `POS`, `REF` and `ALT`. I have 2 VCF files, and I wish to exclude all entries from `VCF1` where it matches those 4 values in `VCF2`. Most of the above-mentioned tools work only on `CHROM` and `POS`, even if they accept a VCF file to filter by. There is no way they compare an exact match to `ALT` allele, even if both VCF inputs were processed to split all multi-allelic sites.

The closest I can get to is by using `bcftools annotate`, and by copying over an `INFO` attribute (say, `INFO/AC`) with a new name (say, `INFO/DUMMY_AC`) so I can filter by that new name. The manual on `bcftools annotate` states:

> When REF and ALT are present, only matching VCF records will be annotated.

which works for me when my `ALT` alleles are split, but does not help me filter, only mark them.

Is there any subset tool that will help me compare by custom attributes or do I have to write my own script for it?

Thank you!

--
Ram

Answer: I've written http://lindenb.github.io/jvarkit/VcfIn.html with the option:

-A, --allalt
ALL user ALT must be found in VCF-database ALT

however, i've not much used this tool.
Question: Potentially a very stupid question, but it wouldn't be my first ever :)

So looking at the sample VCF files from this Biostars (```subset_hg19.vcf```), I see a bunch of lines with ```0|0``` in them.

While I realize this means it has identified that position as matching ```REF | REF``` - why would reference matching regions be in the variants file at all? I thought only variants (positions that DID NOT match the REF in one or both alleles) showed up.

When I run ```bcftools call``` I DONT get any ```0/0```or ```0|0```. Does one have to specifically set the ```call```to identify the matching references - and why would one?

Thanks!

19 416254 rs192385198 T C 100 PASS AC=0;AF=0.000599042;AN=12;NS=2504;DP=18252;EAS_AF=0;AMR_AF=0.0043;AFR_AF=0;EUR_AF=0;SAS_AF=0;AA=T|||;VT=SNP GT 0|0 0|0 0|0 0|0 0|0 0|0
19 416335 rs545360745 G C 100 PASS AC=0;AF=0.000199681;AN=12;NS=2504;DP=18554;EAS_AF=0;AMR_AF=0;AFR_AF=0.0008;EUR_AF=0;SAS_AF=0;AA=G|||;VT=SNP GT 0|0 0|0 0|0 0|0 0|0 0|0
19 416389 rs564036507 G C 100 PASS AC=0;AF=0.000199681;AN=12;NS=2504;DP=17931;EAS_AF=0;AMR_AF=0;AFR_AF=0.0008;EUR_AF=0;SAS_AF=0;AA=G|||;VT=SNP GT 0|0 0|0 0|0 0|0 0|0 0|0
19 416406 rs185752424 G A 100 PASS AC=0;AF=0.000199681;AN=12;NS=2504;DP=17410;EAS_AF=0.001;AMR_AF=0;AFR_AF=0;EUR_AF=0;SAS_AF=0;AA=G|||;VT=SNP GT 0|0 0|0 0|0 0|0 0|0 0|0
19 416449 rs929834 C T 100 PASS AC=0;AF=0.0924521;AN=12;NS=2504;DP=16305;EAS_AF=0;AMR_AF=0.0216;AFR_AF=0.3336;EUR_AF=0.002;SAS_AF=0.0051;AA=C|||;VT=SNP GT 0|0 0|0 0|0 0|0 0|0
0|0

Answer: You would have this if you create a subset of a bigger vcf.

Say that your vcf contains 100 individuals, for which one has a 1/0 call. The rest of the samples will have 0/0.

It looks like subset_hg19.vcf is a, hm, *subset*, of such a bigger vcf ;-)
Question: Usually when I use alignment with bwa, I follow this script:

bwa aln reference.fa sample.fastq.gz > sample.sai

But this time I have manifest.txt and I want to target specific regions with this and I want to speed the alignment process with specific output. If It is possible in BWA, I would like to get your answers. Any help is much appreciated.

Thanks in advance.

Answer: you need to map everything first, convert this "manifest.txt" file as a bed file and extract the regions. https://www.biostars.org/p/49306/
Question: Good afternoon,

I'm using GATK 4.1.7.0 and I'm trying to apply the Haplotypecaller function to a bam file. I was following mainly this post (https://gatk.broadinstitute.org/hc/en-us/articles/360039568932--How-to-Map-and-clean-up-short-read-sequence-data-efficiently) and from that bam file I'm trying to create the GVCF file. Here is my command:

java -jar ~/softwares/GATKK/gatk/gatk-package-4.1.7.0-local.jar HaplotypeCaller --reference Pmuralis_1.0.fa --input mergeandaligned.bam --output mergeandaligned.g.vcf.gz -A StrandBiasBySample -ERC GVCF

I was previously doing the command without the StrandBiasBySample but i saw here (https://gatkforums.broadinstitute.org/gatk/discussion/6813/several-annotations-not-working-in-gatk-haplotype-caller) that was recommended, but I still have the same problem

When I obtain the output, in the log file I have these warnings:

14:54:38.207 WARN StrandBiasBySample - Annotation will not be calculated, genotype is not called or alleleLikelihoodMap is null
14:54:38.658 WARN StrandBiasBySample - Annotation will not be calculated, genotype is not called or alleleLikelihoodMap is null
14:54:38.659 WARN StrandBiasBySample - Annotation will not be calculated, genotype is not called or alleleLikelihoodMap is null
14:54:39.609 WARN StrandBiasBySample - Annotation will not be calculated, genotype is not called or alleleLikelihoodMap is null
14:54:39.610 WARN StrandBiasBySample - Annotation will not be calculated, genotype is not called or alleleLikelihoodMap is null
14:54:41.794 WARN StrandBiasBySample - Annotation will not be calculated, genotype is not called or alleleLikelihoodMap is null
14:54:42.038 WARN StrandBiasBySample - Annotation will not be calculated, genotype is not called or alleleLikelihoodMap is null
14:54:42.038 WARN StrandBiasBySample - Annotation will not be calculated, genotype is not called or alleleLikelihoodMap is null
14:54:43.355 INFO ProgressMeter - Podmur_chrom0001:897923 0.2 3250 19476.6

Later, when I explore the obtained file, I see that for most of the observations, the ALT file is missing, and is as , as you can see here:

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	WA01_1
Podmur_chrom0001	1	.	T	.	.	END=17	GT:DP:GQ:MIN_DP:PL	0/0:0:0:0:0,0,0
Podmur_chrom0001	18	.	T	.	.	END=47	GT:DP:GQ:MIN_DP:PL	0/0:1:3:1:0,3,37
Podmur_chrom0001	48	.	G	.	.	END=48	GT:DP:GQ:MIN_DP:PL	0/0:1:0:1:0,0,0
P

I did not do the next step because I want to solve this first. Thank you very much!

Answer: > Later, when I explore the obtained file, I see that for most of the observations, the ALT file is missing, and is as , as you can see here:

because you are calling in **GVCF** mode (` -ERC GVCF`) : https://gatkforums.broadinstitute.org/gatk/discussion/4017/what-is-a-gvcf-and-how-is-it-different-from-a-regular-vcf
Question: Hi friends

I want to use DESeq2 to normalize the raw count data to do PCA.
I dont have colData.

What code should I use? because in the DESeq2 workflow we need colData and design.
thanks

Answer: You can Google the size normalization algorithm and do it yourself.

Or, make fake coldata, and a design of 1, and then you can make a valid DESeq2 object.
Question: I am getting the following error with fastx_quality_stats with version 0.0.13.2. May be this has to do with the changing Illumina scores. Now they are Phred33 I believe. Has anyone seen this issue

fastx_quality_stats -i test.fastq
fastx_quality_stats: Invalid quality score value (char ')' ord 41 quality value -23) on line 4

Thanks!
-Abhi

more test.fastq
@HISEQ08:235:D142MACXX:7:1101:1474:2097 1:N:0:CGATGT
CCCTCTCCTTCTCCCAGGTAGAACGAGCCTCCCTTTCCCTAGTCAGGGAGTCTTCAAGGAGTTTCTCCTTTGCTGCTGACTTGTCTTTGAATTGCTTGTGGCGTTTTGGAACTTCCTCTTCATTAAGAATCTAGGGTTGTG
CAATTTTCT
+
??<DDD=:C>:CC@A@EEE@C:4:F)@?EDDDD3?B??B9?B?9BD3@@@BC>@CEIIEA=C@?=??. ########################################################################
#########

Answer: A quick google returns: http://seqanswers.com/forums/archive/index.php/t-7596.html

ou need to add the -Q33 parameter to tell it that you're using
Illumina encoded quality scores, not Sanger encoding.