cannin/biostars_qa huggingface dataset:
Question: Hi

I'm having a problem with a whole-exome-sequenced dataset consisting of about 400 human subjects, 200 cases with a certain disease, and 200 controls without. The dataset has been through a rigorous quality control (standardised QC in plink with HWE, IBD, missingness, sex-check ++, along with HapMap population stratification and Eigenstrat/PCA-analysis). I´m using plink to do a basic association-analysis for all variants between cases and controls, and while the resulting QQ-plot for the common (MAF > 0.01) variants is OK, the plot for the rare (MAF < 0.01) variants is less so. Below are the three QQ-plots for all, common and rare variants along with lambda-values:

- QQ-plot all variants, lambda 1.83
http://postimg.org/image/oyr53wchr/
- QQ-plot common variants, lambda 1.04
http://postimg.org/image/6lqjtc20v/
- QQ-plot rare variants, lambda 2.43
http://postimg.org/image/ic4haputb/

The main problem seems to be the positive deviation (observed > expected) of the rare variants in the first part of the plot, causing the lambda to be very big, both for the QQ-plot for rare and all variants. I am wondering what could be the cause of this behaviour for the rare variants, and also what the implications this has for the prospects of doing analysis on rare and common variants together.

I would be grateful if anybody has any experience in these matters and could provide some input.

Many thanks.

Answer: QC has likely nothing to do with your problem. Here are two things you need to think about:

1. you're using a 2df test (genotypic?) for rare variants even though some of your cell counts are likely to be 0.
2. You're using a qq-plot designed to assess the distribution of a continuous variable when yours only has a finite number of possible values due to the small cell counts, making the qq-plot uninterpretable.

Don't over think the qq-plot in that case.
Question: Hello, I have a question for you guys. The perl script rfam_scan.pl uses the program cmsearch to look for a CM model. Wouldn't be better to use cmscan instead of cmsearch? According to the Infernal manual, we should use cmscan if we have a set of sequences and we wanna know the corresponding or closest CM. Don't you think so?

Answer: Hi there, Rfam project leader here (for verification you can check our Twitter feed; we've linked to this question)

Yes, in future we'll be providing cmscan instead of rfam_scan.pl. Historically, rfam_scan.pl needed to include the BLAST filtering step, which reduced the search space the CM needed to cover, making the problem computationally tractable. However, INFERNAL 1.1 now uses HMM-based filters as part of the cmsearches, meaning the blast filters are redundant.
Question: I haven't done Illumina de-multiplxing before as the lab I am used to work with has done this by default. I'm now getting samples from a lab that says that it doesn't do de-multiplexing. It's a new lab which doesn't have experience with NGS so I am not quite sure they now what they are doing. So I got the following Illumina run (this is just a screenshot showing some of the files):

![screenshot of files][1]

---

The SampleSheet in this folder has the following content (again a screenshot):

![samplesheet screenshot][2]

To me it looks like the data has been already de-multiplexed as the fastq files have the same names as the samples. What do you think? How can I test whether this is, indeed, the case? Also, I see some files (their names start with "Undetermined") which are not assigned to any sample. What shall I do with them?

[1]: http://i.imgur.com/tBEef65.png
[2]: http://i.imgur.com/O8HSVLP.png

Answer: These are indeed demultiplexed. Providing demultiplexed files is such a simple thing for a sequencing facility to do that I'm always amazed when one is too lazy to do so. The "Undetermined" files can be ignored. I never provide those to our users. I should note, though, that I routinely check how many reads are in that file, simply because if those files are large then there was an issue with the library or sequencing.
Question: It might be a very stupid question for many of you but, since it's my first variant calling, I didn't figure it out yet.

I have **mpileup**'ped two bam files from two samples, then I filtered the results with **vcfutils.pl** and called the genotypes with **bcftools** **call**. Now I have a **VCF** file containing what I want, but I managed to analyze differences only between the two samples and the reference (which is an assembly coming from a line that is *different* *from both* *samples*).

I have the variants between my two samples and the assembly, but what if I want to detect the differences **between** the two samples? I did it with awk / sed / cut and other command line tools but is there maybe a better and more straight-forward way to do that? Perhaps using bcftools? Up to know, I didn't find it.

Any suggestion appreaciated!

EDIT:
I did it with bcftools gtcheck as well, that works fine. I am asking if there are **other** reliable tools to test!

EDIT 2:
As this post has many views now, I guess it will be useful for everyone to know that I used **bcftools isec** and it worked brilliantly.

Answer: If you are looking for a good way of separating the variants into those that are shared between both samples and those that are unique to each sample, I would recommend using `vcfeval` from [RTG Tools][1] to do the comparison, as that will correctly deal with representation differences between the two samples (this may not be such an issue if both samples were called jointly or you are only concerned with simple SNPs). For example:

rtg vcfeval -t your_reference -b sample1.vcf.gz -c sample2.vcf.gz --squash-ploidy --sample sample1,sample2 -o sample-differences

If both of your samples are in a single VCF, just use the same VCF for both -b and -c. The `--squash-ploidy` option will ignore differences in zygosity in the genotypes and match based on shared alleles. In the output directory you will have variants unique to sample1 in `sample-differences/fn.vcf.gz` and those unique to sample2 in `sample-differences/fp.vcf.gz`, and you can get a summary of the composition of those variants using `rtg vcfstats`.

(disclaimer: I work for RTG, but vcfeval is free and awesome)

[1]: https://github.com/RealTimeGenomics/rtg-tools
Question: Hey,

I have a bunch of SNP-matrices as aligned nucleotide sequences and have no issues in making all kinds of trees out of the alignment. However, I can not find any tool that can plot the actual nucleotide changes (character state changes, mutations) on the branches.
I am dreaming of a tool that can take my .fasta (or .phy) alignment, extract the different states at all positions and plot them on a tree in a "[char_nr]: [change from]>[change to]" way (eg: 12: A>C) or similar to that.

I have tried with the most common classic phylo packages (MEGA7, PAUP4, Mesquite) and also with ETE 3 (python), but all without success (I can get a list of changes and the changes for a single character in PAUP and MEGA, but I am not able to get the information on one tree).
MacClade should be able to do that, but I have no access to a PowerPC-Mac which could run the old code. And WinClada, besides that I don't like spending $50 for trying something, has their servers down at the moment, so no access.

Does anybody have experience with that and can help me out? Thank you so much in advance!

Best

Raphi

Answer: given a tree and its associated alignment, [this script][1] uses ETE3 to plot column changes for each internal branch. Not sure if this is exactly what you need, but it should be easy to add custom changes.
![enter image description here][2]

[1]: https://gist.github.com/jhcepas/868b1057d9c0b72a57b809a67ad6e4cc
[2]: https://cloud.githubusercontent.com/assets/163658/16267090/534671a6-3888-11e6-8eea-24b5c74988d6.png
Question: Hi,

I am trying to filter my GATK 4.0.3 - HaplotypeCaller generated multi-sample VCF for allele depth (AD) annotation at sample genotype-level (so available in "FORMAT" fields of each sample).

I think prior to GATK 4, this annotation was available as "Allele Balance" (AB) ratios (generated by AlleleBalanceBySample), but it is not available anymore in GATK 4. So I tried to filter genotypes based on AD field, that is exactly the same thing but indicated in "X,Y" format, so in an array format of integers. This array format makes it difficult to filter based on depth of alternative allele divided by depth of all alleles at a specific site.

What could be a solution to this? If I could turn this array into a ratio, I could easily filter genotypes using VariantFiltration or other tools such as vcflib/vcffilter. I also tried the below code (following https://gatkforums.broadinstitute.org/gatk/discussion/1255/what-are-jexl-expressions-and-how-can-i-use-them-with-the-gatk):

gatk VariantFiltration -R $ref -V $vcf -O $output --genotype-filter-expression 'vc.getGenotype("Sample1").getAD().1 / vc.getGenotype("Sample1").getAD().0 > 0.33' --set-filtered-genotype-to-no-call --genotype-filter-name 'ABfilter'

This worked, but strangely it filters the variant for all samples if only one of the sample have allele depths that are not in balance (defined by the filter). If it worked only for Sample1, I was planning to write a quick loop for all the samples for instance. I tried the same with GATK 3.8, but still it filters whole variant for all the samples if it is filtered in just one sample.

Any suggestions?

PS: I've asked the same question here also in GATK forum: https://gatkforums.broadinstitute.org/gatk/discussion/12140/allele-depth-ad-allele-balance-ab-filtering-in-gatk-4

Answer: using vcffilterjdk: http://lindenb.github.io/jvarkit/VcfFilterJdk.html

the following expression keep the variants if *all* the genotypes with AD have AD[1]/AD[0]>0.33 .

java -jar dist/vcffilterjdk.jar -e 'return variant.getGenotypes().stream().filter(G->G.hasAD()).map(G->G.getAD()).allMatch(A->A.length>1 && A[0]>0 && A[1]/(double)A[0]> 0.33);' in.vcf
Question: Hi,

Does anyone know how I can modify the command line parameters in needleall to output query coverage? So far my code below only gives me values for Identity, similarity, gaps, and score.

needleall -asequence proteins_f1.fa -bsequence proteins_f2.fa -outfile needle_results -aformat pair -gapopen 10 -gapextend 0.5

Thanks very much.

Answer: A global alignment does not give a coverage. "It will always try a coverage of 100%", a global alignment will try to align the **full** sequence in the best possible way.
Question: Good evening,
I working on whole genome seq. data and I am struggling a bit with plink.
After having merged all my VCF file into a plink .bed file, I need to exclude the SNPs which are not in HW equilibrium.
By merging many whole genome seq VCF files my plink file is full of NA, so I need to encode my SNPs as 0/1/2 (to set the missing as homozygous for the reference allele) before verifying the HW condition.
After that I need to recode my file as traw file (as an output from plink --recode A-transpose, which has the individuals in columns and a line for each SNP).
I would need something like recode12 --fill-missing-a2, I guess, but recode12 output just a ped (or a tped) file, which is not what I need.
Does anyone know a possible solution?
In few words I need to encode my plink .bed file as 0/1/2 (2 stading for homozygosity for the minor allele, and 0 to fill all the NA), in order to check the HW eq. and then convert it to a matrix with an individual fro each column and a SNP for each line).
Thanks

Answer: the way I do it is to first genotype the vcf files I have with a common set of SNPs (in other words genotype against a list of SNP positions, so that every individual genotyped is going to have the same number of markers).

bedtools intersect -a my.vcf -b bed_present_snps -wa -header > out.vcf

the annotate the file:

cat my_out.vcf | vcf-annotate -a vcf_annotation_file_sorted.gz -c CHROM,FROM,TO,ID,-,-,- > annotated.vcf

and move it to plink format:

./plink --double-id --vcf my.vcf --recode --out my_plink

from there it's just converting to binary and merging with other files / filtering:

./plink --file my_plink --make-bed --out plink_binary

./plink --file file_with_all_snps --bmerge my_plink_binary --make-bed --out plink_merged

now the missing markers should already be set as 0 0, and you can recode however you want to suit your liking
Question: Hi, I am looking to import .nexus file in R to analyze it but am facing errors in the import step itself. I tried the following commands and obtained the following error messages.

```
> file<-read.nexus("hominoid.nex")
Error in start:end : NA/NaN argument
```

```
>tree<-read.tree("hominoid_apriori.tre")
Error in if (tp[3] != "") obj$node.label <- tp[3] :
missing value where TRUE/FALSE needed
```

```
> file2<-read.nexus("testmodel.nex")
Error in X[sb] <- gsub(".*\\]", "", X[sb]) :
NAs are not allowed in subscripted assignments
In addition: Warning messages:
1: In grep("\\[", X) : input string 26 is invalid in this locale
2: In grep("\\[", X) : input string 231 is invalid in this locale
3: In grep("\\]", X) : input string 26 is invalid in this locale
4: In grep("\\]", X) : input string 231 is invalid in this locale
5: In LEFT == RIGHT :
longer object length is not a multiple of shorter object length
```

The file link,

https://www.dropbox.com/sh/kqo8bwi9s7ispcx/AADDEsRuDfLCpKbZeGgU81nha?dl=0

I will be glad to provide necessary details, kindly advice.

PS: I am able to read the file using `read.csv` command but not able to perform any operations with it.

BR, M

Answer: In ape, `read.nexus()` reads in a treefile formatted in NEXUS format. `read.tree()` reads a treefile formatted in Newick (or parenthetical) format.

You want `read.nexus.data(file=hominoid.nex)` for the data and `read.nexus("hominoid_apriori.tre")`.

`testmodel.nex` is just a PAUP* block and not a properly formatted NEXUS data file. This is appended to a NEXUS data file and then executed. If you want to do model selection, may I suggest Minin et al.'s decision-theroetic approach, [DT-ModSel][1]? Append the block, execute the NEXUS file, then run the perl script.

Please respond if things aren't working.

[1]: http://www.webpages.uidaho.edu/~jacks/DTModSel.html
Question: Hello,

I would like to know how the string of the field "LOCUS" on a Genbank (and other similar) entry is formed, and how to read its information. For example: what does NM stand for? What about the XR or XM entries? How do the number following the letters is obtained?

Answer: The meaning of the LOCUS field is described [here](http://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html#LocusB)

The NM and XM entries, etc., are RefSeq identifiers rather than GenBank identifiers. RefSeq is effectively a well-curated subset of GenBank.

[This page](http://www.ncbi.nlm.nih.gov/RefSeq/key.html) explains the meanings of the various prefixes.
Question: Hi All -

I'm using the R package PopGenome to analyze WGS data from a VCF file. I read in my VCF and GFF files by chromosome, as described [in this post][1]. Now I'd like to verify synonymous and non-synonymous substitutions using the `set.synnonsyn()` function. The `set.synnonsyn()` function has a `ref.chr = ` option, but I cannot find any details on what format the input chromosome should take when there are in fact multiple chromosomes in play. Should I feed the program a multifasta? Should I point it to a directory that contains individual chromosome fastas, similar to the way we import multiple VCF and GFF files? I've tried this every way I can think of but no luck.

Any help is much appreciated -

[1]: https://www.biostars.org/p/121428/

Answer: Hi,

test <- set.synnonsyn(test, ref.chr=c("chr1.fas","chr2.fas", ... ))

should work fine.

Best
Question: I am not interested in read mapping *per se*.

I am interested in reporting **ALL** possible end-to-end alignments of a read against a reference that have <= k mismatches (e.g. k = 5). I do not care about the "best" alignment.

Standard DNA & RNA read aligners (bwa, bowtie, STAR, etc.) that try to find the "best" alignment for a read place limits on number of mapping locations for short seeds, or something similar.

Does any read aligner allow FULL sensitivity?

For example, I can use MUMmer3 to report **ALL** exact matches of a read in a genome (of a given substring length). I am looking for something similar - i.e. full sensitivity - but allowing a limited number of mismatches.

Anybody know of a program that does this? I have been looking but am not convinced of the full sensitivity of the aligners.

Answer: There are several "all" mappers, [segemehl][1] is one I see recommended frequently. Also, see [this][2] (long but even so incomplete) list, search for "A" under "alignment reported" for all-mappers. Finally, if you are wary of results, create a small test data set and benchmark the tools.

[1]: http://www.bioinf.uni-leipzig.de/Software/segemehl/
[2]: http://www.ebi.ac.uk/~nf/hts_mappers/
Question: Hi,

I was wondering if there was a way to sort a fasta file so that the chromosomes are in numeric order. Currently, the fasta file is computationally sorted where the chromosomes are order like this

>chr1
>chr10
>chr11
....
>chr2
>chr20

I was wondering if there was a way to sort the fasta file so the chromosomes are sorted numerically like this

>chr1
>chr2
>chr3
...
>chr10
>chr11

Answer: Try:

bioawk -c fastx '{print}' in.fa | sort -k1,1V | awk '{print ">"$1;print $2}'

or any of:

https://www.biostars.org/p/18929/
Question: Hi,

I searched the previous posts and the nanopore community forum but couldn't find a definitive answer so here I am looking for your advice. I know working with long-reads RNA-seq data is still complicated because it's hard to find a general consensus on how to analyse them but I'd be interested in having your opinion.

I'm interested in measuring the number of genes detected in various RNA-seq libraries that I have generated with various nanopore library kits (DNA ligation of RT-PCR products, direct-cDNA and direct-RNA).

So far I've been using **HTseq** with the mode `intersection-nonempty` but it seems that now the `union` mode is the most recommended one. And i've also seen posts stating that **featureCounts** is much better and there's even an option for `longreads`.

I tried to compare the output of the three 'possibilities' but I cannot say it really helped me in the end:

- HTseq + intersection-nonempty mode: 238K assigned + 0 ambiguous + 18K no_features + 24K too low quality

- HTseq + union mode: 210K assigned + 44K ambiguous + 2K no_features + 24K too low quality

- featureCounts + long-reads mode: 219K assigned + 59K ambiguous + 2K no_features + 0 too low quality

Based on those observations, would you be able to recommend me to use one versus the two others ? And if yes, why ?
I'm trying to understand what makes one more suited to my need and/or generally more suited to the analyse of long reads RNA-seq data.

Thanks a lot.

Answer: Having thinking a little more about the problem, I think you may try Salmon with its alignment-based mode of quantification. The expectation-maximization algorithm will probably do a better job at sorting ambiguous reads than HTSeq / featureCounts. Be sure to read the manual: https://salmon.readthedocs.io/en/latest/salmon.html
Question: Hey all :)

The SAM spec is pretty lenient on what is and is not allowed in a QNAME. The only really relevant part is:

> QNAME: Query template NAME. Reads/segments having identical QNAME are
> regarded to come from the same template. A QNAME ‘*’ indicates the
> information is unavailable. In a SAM file, a read may occupy multiple
> alignment lines, when its alignment is chimeric or when multiple
> mappings are given.

So each ***template*** gets a unique QNAME, and the same QNAME can appear multiple times in the file (for multiple sequenced reads, secondary alignments, etc).

The problem i'm facing right now is that there's obviously a ***lot*** more to QNAMEs than just that. First off, I thought a template was not a read, but the whole fragment as it entered the sequencer - so paired-end reads are two reads from the same template, so they should have the same QNAME. In other words, putting /1 and /2 or anything else at the ends of the QNAME to denote the read pair should go against the standard?

The QNAMEs often also encodes the position in the flowcell that the template came from, which some programs use to detect optical duplicates. However, I suspect that duplicate detection is exactly the same without it - just those duplicates are marked as PCR duplicates rather than Optical (which is information usually thrown away anyway).

So my question is, if I was to rename all the QNAMEs in a BAM to something unique to the template, but without the flowcell info nor the /1 /2 mate info, would that cause any actual problems downstream? Are there any tools that just wont work if I just rigorously follow the standards here? If so, what would a good 'fake' QNAME look like?

Answer: Short answer: read names don't matter except for sorting (by read name, not coordinate), pairing (e.g., counting with featureCounts or htseq-count), and marking optical duplicates. Of course optical duplicates will get marked as duplicates regardless, so who really cares about that.

BTW, typically aligners strip /1 and such off, though not always. One should generally not rely on qnames for anything unless absolutely needed.
Question: Hello,

I have multiple fasta sequences that are like this:

>2p__scaffold_2__5799__6580__-__778568__0.00__0.00
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGCAGCCGGA
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAA
>2p__scaffold_2__5799__6580__+__778569__0.00__0.00
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGC
>1p__scaffold_2__11235__11438__-__830827__0.00__0.00
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGCAGCCGGA
GCTTCAATCCAGGGGATCGAGGAGATCCAAAGCAGCAGAAGCGGCTCGACGATGGTGAGG
ATTCGGGATCGGATTCAGCGCTCGTCGGGACTGG
>1p__scaffold_2__33129__34129__+__811706__0.00__0.00
GCTGGCGACGGATCTA

And I want to keep just the "> + ID" (numbers after `__+/-__` and before __0.00_0.00)

So I expect an output like this:

>778568
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGCAGCCGGA
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAA
>778569
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGC

I searched for it and tried this:

sed 's@.*__-__@@' input.fa > output.fa

That removed `__-__` and everything before it, including the ">" that I wanted to keep.

I also tried this to remove everything between ">" and `__-__`

sed -e 's/\>//' -e 's/\__-__.*//' input.fa > output.fa

But this removed everything after `__-__`

And this, that removed __0.00_0.00

sed 's/__0.00.*$//' input.fa > output.fa

Thank you for your help.

Answer: Try the following regular expression:

$ awk '{ if ($0 ~ /^>/) { match($1, /[+|-]__[0-9]+__/, m); print ">"substr(m[0], 4, length(m[0]) - 5); } else { print $0; } }' input.fa > output.fa

Then:

$ less output.fa
>778568
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGCAGCCGGA
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAA
>778569
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGC
>830827
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGCAGCCGGA
GCTTCAATCCAGGGGATCGAGGAGATCCAAAGCAGCAGAAGCGGCTCGACGATGGTGAGG
ATTCGGGATCGGATTCAGCGCTCGTCGGGACTGG
>811706
GCTGGCGACGGATCTA