cannin/biostars_qa huggingface dataset:
Question: I'm kind of new to this space-- a friend of mine says he uses SNPeff for all his exome annotations, and he doesn't know of any other popular tools for this purpose.

I'm annotating some human exomes and I am curious about what else is out there. A search gave me a lot of answers, but I don't know which are popular in the community. Are there gaps the SNPeff leaves that other effect predictors fill? Thank you so much for reading my post!

Answer: The tools I hear used most frequently are SnpEff, VEP, and Annovar. This [paper][1] (Table 1) shows a comparison of the three tools.

SnpEff tends to be robust and I personally use it the most. Remarkably, SnpEff can effectively annotate even structural variants and long indels, in addition to traditional smaller variants. I've used Annovar once or twice but strange bugs crop up here and there; however the developer of it maintains it well and offers a lot of [documentation][2]. VEP seems quite popular, but I personally have the least experience with this one.

[1]: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0974-4
[2]: http://annovar.openbioinformatics.org/en/latest/
Question: How are variant callers able to compute which copy of the chrs in case of a diploid species has the heterozygous detected variants? What information do they use, just paired end reads?

I could find some statistical data in GATK webpage, but I would like to understand if there is other information used, the rationale behind it the accuracy it would have and the factors that affect this process.

Thank you very much

Answer: Usually the caller has no idea which chromosome homolog a variant is on. It can just see variants that are in the same read or read pair (unlikely for short reads) or it can try to infer which variants are on the same chromosome homolog (phased) using read-backed phasing (as part of the read assembly performed by the haplotype caller).

These in silico methods are spotty at best. Most people who need phasing just use a long-read technology, or they sequence the parents.

https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_phasing_ReadBackedPhasing.php

https://www.illumina.com/techniques/sequencing/dna-sequencing/whole-genome-sequencing/phased-sequencing.html
Question: I have a dataset of 32 samples quantified using `salmon`. I used `tximport` to import the data:

txi=tximport(c(a vector of file names),type="salmon",tx2gene=tx2gene)

I have created the colData and used `DESeqDataSetFromTximport` to construct the Deseq2 data set:

dds <- DESeqDataSetFromTximport(txi = txi,colData = col_data,design = ~ Form)

In this case, all 32 samples are used, now I only want to use a subset of txi by choosing specific samples. Is there a way to subset my current txi rather than constructing a new txi using only the samples that I want?

Answer: The DESeq2 object is basically a SummarizedExperiment so you can subset this prior to running `DESeq` using standard operations such as `dds_use <- dds[,c(columns_to_keep]` where `columns_to_keep` is a numeric vector with the columns (samples) you want to use.
Question: I'm trying to produce a simple window => coverage table from a bam file. I produced (non-sliding) windows with `bedtools makewindows` and then ran `bedtools coverage -abam alignment.bam -b reference.windows > coverage`. The problem is that I cannot make sense of the output I got. Here's a sample:

```
I 0 215 SRR1297046.176411/2 23 + 0 215 0,0,0 1 215, 0, 1 215 215 1.0000000
I 0 216 SRR1297046.2210768/1 23 + 0 216 0,0,0 1 216, 0, 1 216 216 1.0000000
I 1 222 SRR1297046.178368/1 24 + 1 222 0,0,0 1 221, 0, 1 221 221 1.0000000
```

The first field is the chromosome name and the fourth one is the read's, but what are the others? This output is inconsistent with the ones described [here][1]. How can I convert this output into a "window => coverage" table?

I'd appreciate any help.

[1]: http://bedtools.readthedocs.org/en/latest/content/tools/coverage.html?highlight=coverage

Answer: If you want the coverage of each entry in reference.windows then bedtools coverage -a reference.windows -b alignment.bam > coverage is what you want.

BTW, this will probably double count overlapping paired-end reads. You could alternatively use samtools depth with some post-processing and not have to worry about that.
Question: How do I get the start and end positions of the Hemoglobin Z transcript ENSG00000130656 using a Python REST API? Is it possible? For example I want to get the start and end positions of these exons: http://www.ensembl.org/Homo_sapiens/Transcript/Exons?db=core;g=ENSG00000130656;r=16:152687-154503;t=ENST00000252951

Answer: I came up with the answer. I can simply access this url: http://rest.ensembl.org/overlap/id/ENSG00000130656.json?feature=exon

The documentation is quite comprehensive, once it is located: http://rest.ensembl.org/documentation/info/overlap_id
Question: Are all [GO terms][1] manually curated or are they computationally curated or mixture of both?

How about GO from [MSigDB][2] ?

[1]: http://www.geneontology.org/
[2]: http://software.broadinstitute.org/gsea/msigdb/

Answer: As far as I know, the ontology itself is entirely curated. On the other hand the gene annotations are a mixture of curation and automatic annotations. This is tracked with [evidence codes][1].

[1]: http://www.geneontology.org/page/guide-go-evidence-codes
Question: Hi all,

I have annotated my genome with the classical combo MAKER + Blast2GO. I am writing now a paper about this genome and this protein coding gene annotation.

What is it usually do in term of protein coding gene annotation upload ? Do people upload it in their github or specific database ?

Thanks for your help.

Answer: For the European entry to INSDC use ENA.
For the US door to INSDC use NCBI.

In both cases you must convert your annotation.
For submission to ENA you need to use [EMBLmyGFF3][1].
For submission to NCBI you need to use [GAG][2].

[1]: https://github.com/NBISweden/EMBLmyGFF3
[2]: https://github.com/genomeannotation/GAG
Question: Hi,

I have a file with five columns, from which I want the mean of the last four columns where the first rows of the first column are the same. For example, my file looks like this -

A 1 1 1 1
A 1 3 1 1
A 1 1 2 3
B 5 7 2 4
C 2 1 5 1
C 2 2 3 6

The desired output is -

A 1 1.7 1.3 1.7
B 5 7 2 4
C 2 1.5 4 3.5

Any help would be appreciated.

Thanks.

Answer: import pandas as pd
df = pd.read_csv("tmp.txt", sep="\t", header=None)
df
0	1	2	3	4
0	A	1	1	1	1
1	A	1	3	1	1
2	A	1	1	2	3
3	B	5	7	2	4
4	C	2	1	5	1
5	C	2	2	3	6

df.groupby(0).mean()

1	2	3	4
0
A	1.0	1.666667	1.333333	1.666667
B	5.0	7.000000	2.000000	4.000000
C	2.0	1.500000	4.000000	3.500000

df.groupby(0).mean().to_csv("tmp_mean.txt", sep="\t", header = None)
Question: I'm looking for a Python package that has all the capabilities of **METATOOL 5.0** for doing **Elementary Mode Analysis**. I found **ScrumPy** and **PyNetMet**. I'm having a little trouble installing all the dependencies for ScrumPy and before I proceed, I was wondering if it was worth it for doing what I want to do. Any help would be greatly appreciated :)

If you're not familiar with Elementary Mode Analysis, "These are the smallest sub-networks that allow a metabolic reconstruction network to function in steady state.[[28]][1][[29]][2] According to Stelling (2002),[[29]][2] elementary modes can be used to understand cellular objectives for the overall metabolic network. Furthermore, elementary mode analysis takes into account [stoichiometrics][3] and [thermodynamics][4] when evaluating whether a particular metabolic route or network is feasible and likely for a set of proteins/enzymes.[[28]][1]" from the wiki http://en.wikipedia.org/wiki/Metabolic_network_modelling#Elementary_mode_analysis.

[1]: http://en.wikipedia.org/wiki/Metabolic_network_modelling#cite_note-Schuster_2000-28
[2]: http://en.wikipedia.org/wiki/Metabolic_network_modelling#cite_note-Stelling_2002-29
[3]: http://en.wikipedia.org/wiki/Stoichiometry
[4]: http://en.wikipedia.org/wiki/Thermodynamics

Answer: **PySCeS** is the best option for there are regular updates and more documentation then you would ever need to read.
Question: Hi

I would like to plot RNASeq data that I have downloaded from TCGA in a PCA plot. I have found some great guides on how to plot the actual data in PCA using r in ggplot2 and such but my main question is what format data should I plot?

I currently have raw counts and RSEM data. Should I input raw counts into something like edgeR or deseq2 and filter for expression by cpm first? Should I normalise it? Should I stabilise variance using rlog2? Or convert to TPM and plot that? Argh I'm so confused. Grateful for any advice you can give me :)

Answer: Dear Elizabeth,

In the simplest scenario (4 samples; 4 genes; 1 experimental condition), your metadata object, which you may have to read in from a file, could look like:

ID	Condition
MA	M
MB	M
PA	P
PB	P

...whilst your counts file could look like:

MA	MB	PA	PB
gene1	45	46	25	22
gene2	45	45	45	44
gene3	10	10	9	4
gene4	88	67	34	44

This could then be read into DESeq2 as:

dds <- DESeqDataSetFromMatrix(rawcounts, colData=metadata, design=~Condition)
Question: Hi,
I have a set of files that I’d like to perform a function on, with the goal of applying one or more parameters in that function that include more than one possible state.

For example, I might have two samples, each with their own fasta file: sample_A and sample_B.

I want to perform a blast search for each input fasta file, but I also want to loop through a range of word sizes for every blast process for each sample. Say, three values: 11, 13, 15.

This would mean that for the sample_*.fasta input, I’d generate three blast output files, each one reflecting one of those three word size values.

I am struggling to understand how to structure the snakemake rule for input and output names, because my they don’t share the same wildcards - there is an extra name from the blast parameter output that isn’t part of the input name.

Thanks for advice on how to include a parameter name in a snakemake rule into the output name!

Answer: Maybe this?

```
samples = ['sample_a', 'sample_b', 'sample_c']
word_sizes = [11, 13, 15]

rule all:
input:
expand('blast/{sample}.{word_size}.out', sample= samples, word_size= word_sizes),

rule blast:
input:
fa= '{sample}.fa',
output:
out= 'blast/{sample}.{word_size}.out',
shell:
r"""
blastn -word_size {wildcards.word_size} -query {input.fa} -out {output.out} ...
"""
```

Note that the `expand()` function will create all combinations of sample and word_size and returns a list of strings. If you want more control on what combinations to have you can use any python code to create such list.

Also, this assumes fasta file are named with the `sample` prefix, if this is not the case you can use a dictionary to map samples to fasta files. In this latter case you may need to use a function as input to the `blast` rule.
Question: Hi,

I have a list of variants called from a individual genome and I'm trying to filter out the important predispositions from it. My approach was to download the `variant_summary.txt.gz` file from ClinVar website, in which most of the variants related to human health are being recorded, so that I can intersect my variants with it.

I loaded the `variant_summary.txt` into R and it says the Dataset has 154358 rows and 25 columns. But when I check with `wc-l` linux the number of records is 198661. I double checked the no of rows by visualizing the data in excel. It had 198,661. My questions are,

1. Why R does not load all the records of my file?
2. Given the fact that I'm still novice to bioinformatics do you think that my approach is feasible in finding predispositions if I fix the R issue?

Thanks you very much.

Answer: You can read in the file like so:

```
var.anno = read.delim("Downloads/variant_summary.txt", header=T)
dim(var.anno)
[1] 198660 25
```

That gave me the correct dimensions, first row contains the header.

Regarding your second question, I think it is reasonable to try to annotate detected variants and their association with phenotypes from available databases. Depending on what you are after, you might also want to consult other variant databases, e.g. dbSNP. dbSNP also links to ClinVar if there is [an entry there][1] and to OMIM.

There are also R packages for this purpose, some here:

- http://cran.r-project.org/web/packages/NCBI2R/index.html
- SNP-related packages in bioconductor http://www.bioconductor.org/packages/release/BiocViews.html#___SNP

[1]: http://www.ncbi.nlm.nih.gov/SNP/snp_ref.cgi?rs=4988235
Question: Greeting all.

I am searching novel somatic mutations using various public software and annotation tools such as annovar.

I have been using my own cancer exome data not public data 'just real patient data' and finally I completed running calling program(MuTect and annovar) and I got to discover some genes that looks disease-like to me.

One of my genes is starting with LOC and its full name is LOC100129697 and it has variation in nonsynonymous region.

I hope that this gene are related to the disease and with this hope in mind, I ran the IGV to analyze this variation more deeply.

However, when I looked up IGV , I was surprised to see that gene name which IGV shows me is totally different from LOC100129697.

How can I interpret this gene LOC100129697?

Is there a gene which have another nick-name?

Answer: From NCBI:

> **Symbols beginning with LOC.** When a published symbol is not available, and orthologs have not yet been determined, Gene will provide a symbol that is constructed as 'LOC' + the GeneID. This is not retained when a replacement symbol has been identified, although queries by the LOC term are still supported. In other words, a record with the symbol LOC12345 is equivalent to GeneID = 12345. So if the symbol changes, the record can still be retrieved on the web using LOC12345 as a query, or from any file using GeneID = 12345.
Question: Hi,

I would ask if someone could clarify me the concept of phased/unphased genotypes. In the first case we know the haplotype and in the second case we don't know the haplotype ?

Thanks.

Answer: You know the haplotypes in each case. The difference is that for phased heterozygous variants, you know which variants are on the same chromosomal copy. For example, take heterozygous variant 1 (A/C) and variant 2 (T/G). One chromosome copy will have an A at location 1, and the other will have a C. But what allele will be present at location 2 for the chromosome copy with an A at location 1? With unphased variants, you don't know - could be T or G. With phased variants, you know which one it is. Phasing is not necessarily complete; often you can only phase small regions.
Question: Hi all,

As you see in the picture, I have two columns. So, I want to replace the values smaller than 0.500000 from the column A_Freq with the values in the M.F column In the same rows. What is the best idea?

for example: replace 0.312500 In the third row of the column A_Freq with 0.687500 From the same row but in the column M.F.

CHROM_POS A_Freq M.F N_Chr
1 CM009840.1_1096 0.812500 0.187500 16.25000
2 CM009840.1_1177 0.611111 0.388889 12.22222
3 CM009840.1_1276 0.312500 0.687500 6.25000
4 CM009840.1_1295 0.277778 0.722222 5.55556
5 CM009840.1_1471 0.250000 0.750000 5.00000
6 CM009840.1_1518 0.875000 0.125000 17.50000
7 CM009840.1_1527 0.222222 0.777778 4.44444
8 CM009840.1_1533 0.777778 0.222222 15.55556
9 CM009840.1_1630 0.250000 0.750000 5.00000
10 CM009840.1_1639 0.000000 1.000000 0.00000
11 CM009840.1_1711 0.500000 0.500000 10.00000
12 CM009840.1_1972 0.250000 0.750000 5.00000
13 CM009840.1_2030 0.142857 0.857143 2.85714
14 CM009840.1_2101 0.375000 0.625000 7.50000
15 CM009840.1_2690 0.687500 0.312500 13.75000
16 CM009840.1_2849 0.142857 0.857143 2.85714
17 CM009840.1_3013 0.312500 0.687500 6.25000
18 CM009840.1_3042 0.714286 0.285714 14.28572
19 CM009840.1_3062 0.250000 0.750000 5.00000
20 CM009840.1_3128 0.250000 0.750000 5.00000

Best Regard

Modtafa

Answer: library(tidyverse)
new_df <- mutate(df, A_Freq = ifelse(A_Freq >= 0.5, A_Freq, M.F))
Question: Dear Colleagues,

I would like to do some visualization of my network in Cytoscape. My file contains genes and their expression values in normal and tumor samples. The visualization of interest would be the same network colored according to A) expression in normal samples and B) expression in tumor samples.

Honestly, I completely do not know how should be prepared the input file with such annotation. Basic input file to create network in Cytoscape looks like:

Gene A - Gene B
(and some additional annotations).

Some example of such analysis shown in Cytoscape tutorial indicates that such input file should look like this:

GENE COMMON gal1RG gal1RG gal80R

YHR051W COX6 -0.034 -0.034 -0.304

YHR124W NDT80 -0.090 -0.090 -0.348

YKL181W PRS1 -0.167 -0.167 0.112

YGR072W UPF3 0.245 0.245 0.787

(unfortunately, it is not described how the values were obtained).

The above arises in my question: what these values are related to?
If in my file I will put expression values for genes only in first columns, then I won't be able to color genes from second column. Moreover, in Cytoscape I can use only one column as such annotation, so there is no chance to provide two columns with expression values for both nodes (start and end).

I would appreciate any suggestion how to solve such issue. Maybe there is some mathematical solution how to calculate "something" that I could possibly use as representation of expression or maybe there is some other software that would be able to do something like this?

Thanks in advance.

Best regards!

Answer: **INPUT DATA**

You need two input files to color your network according to expression values in normal and cancer samples.

First file defines interactions between genes. For example, Excel file (`interactions.xls`):

Source	Target
Gene1	Gene2
Gene1	Gene3
Gene2	Gene3
Gene2	Gene4
Gene3	Gene4
Gene3	Gene4
Gene4	Gene5

Second file lists gene names with their expression values in normal and cancer samples. For example, Excel file (`expressions.xls`):

Gene_name	Normal	Cancer
Gene1	1.92	1.56
Gene2	0.56	2.02
Gene3	0.34	0.01
Gene4	0.98	0.98
Gene5	0.02	0.87

**CYTOSCAPE**

1. Open Cytoscape.
2. `File` > `Import` > `Network` > `From file`
3. Load the `interactions.xls` file. Define first column as *Source*, and second as *Target*.
4. Load the `expressions.xls` file.
- `File` > `Import` > `Table` > `File`.
- Label first column as *Gene_name*, second as *Normal* and third as *Cancer*.
5. Color your network according to expression values in the *Normal* column.
- `Style` > `Fill color` > `Column`: Normal
- `Mapping Type`: Continuous
- `Current mapping`: choose your colors.

**OUTPUT**

![enter image description here][1]

Pozdrawiam:)

[1]: http://i66.tinypic.com/2dv770k.png
Question: Dear all,

I have a molecule in PyMol which I want to open in ChimeraX. But after exporting the molecule as .pdb, some structural information seems to be lost - regions that in the PyMol session file are embedded into helices are now loops (this regards to several residues only). If I open the .pdb in PyMol, it looks the same as in chimeraX. I have no idea what I keep doing wrong, any hints? The first 14 residues of the original protein (from PDB, corresponding to Histag) were removed with a custom PyMol script which changed residues enumeration so that it properly starts with 1 - that is why the molecule from PyMol is needed. I would be grateful for the help.

best,

maciek

Answer: For anyone who might be interested:
1 - open original .pdb in text editor
2 - remove all atoms (lines starting with ATOMS) corresponding to unwanted residues, save file
3 - go to https://wenmr.science.uu.nl/pdbtools/ and renumber your sequence (saved file) as desired
4 - go back to your original structure file and find all lines starting with HELIX or SHEET and copy them to the beginning of your corrected file (just before the first residue) and apply corrections to the secondary elements ranges (for example, I deleted 14 residues from the beginning, so need to shift back secondary elements ranges by 14)
5 - save the file and check if the structure adheres to what you expect

For me, it worked perfectly.
Question: Hi,

I've a bigwig containing a chip-seq signal for a histone mark (H3K4me1) and I want to produce a bed file containing the peaks (the same way MACS produced a bigwig and a peak bed file)? How can I do that? If it's impossible I can realign the raw data and do the whole peak calling step but I prefer to avoid that.

Thanks

Answer: MACS2 can call peak from bedgraph, you can convert bigwig to bedgraph first by bigwigtobedgraph from [here][1] and then call peak with [MACS2][2].

> bdgpeakcall: Call peaks from bedGraph output.

[1]: http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/
[2]: https://github.com/taoliu/MACS/
Question: Hello,
I have a fasta file with lines of format -

>FBti0019256 type=transposable_element; loc=2L:22300300..22304444; name=invader2{}555; dbxref=FlyBase_Annotation_IDs:TE19256,FlyBase:FBti0019256; MD5=d9259a0e33aad699215e64916bd47a5b; length=4145; release=r6.19; species=Dmel;

I would like to convert these lines into a bed file of format -

chr2L /t 22300300 /t 22304444 /t invader2

Is there a program that can directly perform this conversion or an awk command that can do this easier? Please let me know, thank you for your help.

Answer: use an associative array in awk to store each component.

I'm just too lazy to extract the chrom/start/end, you get the idea.

awk -F ' ' '/^>/ {delete map;for(i=1;i<=NF;i++) {eq=index($i,"=");if(eq==0) continue;key=substr($i,1,eq-1);val=substr($i,eq+1); gsub(/;$/,"",val); map[key]=val;} OFS="\t"; print map["loc"],map["type"],map["name"],map["dbxref"],map["MD5"];}'
Question: Hi all,

I have paired-end 4 ATAC-seq data (2 replicates for 2 samples). I have done aligning using Bowtie2. I did filter MT reads and duplicates using Picard, then performed peak calling on Bam file using MACS2. Also I did differential peak analysis using deeptools and filter them by FDR<0.05 and abs(2foldchange)>2.

After these, I generated density peak heatmaps using deeptools. However, on the figure top on the heatmaps height of peaks are not the same for 4 files although I normalized bam files while converting to bigwig using bamCoverage.

My questions are:
Should the height of those peaks be the same or slight change is acceptable?
If not how can I normalize the data? Should I normalize bam files then do the peak calling again if so which tool you suggest? or diffbind normalization okay?
Also, I am really confused about the coverage file normalization and peak normalization.
Lastly, as written in this post https://www.biostars.org/p/308976/ how can I downsample each sample?

If you could explain these I will appreciate it.

![peak signal heatmap][1]

[1]: https://www.dropbox.com/home?preview=image.png

Answer: I wouldn't expect the peak heights to be identical, some amount of biological variation is normal. The goal of the normalization should instead be to set the background level to roughly similar values between samples.

You do not need to normalize your BAM files, CSAW or DiffBind will take care of that step for you.

If you want to downsample you can either use `samtools view -s` if starting from BAM files or `seqtk` if starting from fastq files. This is generally not needed.