cannin/biostars_qa huggingface dataset:
Question: Hi,

I have VCFs of different samples split by chromosomes. I want to merge all these chromosomes files. For this I am using this command

bcftools merge --force-samples -m all *vcf.gz

The size of all vcf in gz form is 230GB. But the file which is produced after merging is of 1.2TB. And it has merged only two chromosomes till now. About 16 hours has passed when I ran this command.

Is something wrong with the command or it is running normally?

Answer: I believe it is expected that the merged size will be larger than the sum of all individual VCFs. Why? - This is the case because, in the merged file, a record is created for every single variant record across all of your files, and, even if the variant is only present in a single sample (private variant), information about it has to be added for each individual.

Consider the following:

**Sample1**

chr1 1234 A T 1/1
chr1 5678 A G 0/1
chr1 9999 G C 0/1
chr6 1983 T C 0/1

**Sample2**

chr3 1000 A G 0/1
chr6 1983 T C 1/1

**Sample3**

chr1 4444 T C 1/1

**Merged:**

S1 S2 S3
chr1 1234 A T 1/1 ./. ./.
chr1 4444 T C ./. ./. 1/1
chr1 5678 A G 0/1 ./. ./.
chr1 9999 G C 0/1 ./. ./.
chr3 1000 A G ./. 0/1 ./.
chr6 1983 T C 0/1 1/1 ./.

The situation becomes 'worse' the more heterogeneous is your cohort, and also if your samples have many private variants.

Coincidentally, if you save your merged file as BCF, it will likely be less than the total size of your individual VCFs because BCF is binary format. Think of the difference between SAM and BAM - huge difference in size.

If you still have major issues with space and time, then consider removing all private variants prior to merging, or those variants with MAF lower than a certain cut-off.

----------------------------------
----------------

Before merging, by the way, I would consider normalising each and every individual VCF with:

while read gzVCF
do
bcftools norm -m-any "${gzVCF}" | bcftools norm -Ob --check-ref w -f human_g1k_v37.fasta > "${gzVCF}".bcf ;

bcftools index "${gzVCF}".bcf ;
done < gzVCF.list

- 1st pipe, splits multi-allelic calls into separate variant calls
- 2nd pipe, left-aligns indels and issues warnings when the REF base in
your VCF does not match the base in the supplied FASTA reference
genome

The file gzVCF.list contains the path of all of your VCFs / VCF.gz. This will ensure that your merge is conducted faithfully.

Kevin
Question: I have 3 replicates data outputs in separate three .csv files.E

Each files the first column is Gene identification numbers (AGI: my experimental targets). The three files consisted of different experimental outcomes. i.e. first file 250 targets, second ~500, and third ~ 2000.

Each file consisted with other information for the identified targets, e.g.

I want to identify common targets (AGIs) in these three .csv files and make a new working table combining other details of the common AGIs.

How can I identify the common targets (AGIs) using R programming?

Answer: Using ***tidyverse***, try this example

# dummy data
rep1 <- data.frame(AGI = c(1:5,6), value = 1)
rep2 <- data.frame(AGI = c(2:5,6), value = 2)
rep3 <- data.frame(AGI = c(1,5,6), value = 3)

library(tidyverse)

bind_rows(list(rep1, rep2, rep3), .id = "Rep") %>%
spread(key = "Rep", value = "value", sep = "") %>%
filter(complete.cases(.)) %>%
mutate(AverageValue = rowMeans(.[, -1]))

# AGI Rep1 Rep2 Rep3 AverageValue
# 1 5 1 2 3 2
# 2 6 1 2 3 2

---------
**Tip:** We could create the list of dataframes automatically, something like:

myReps <- lapply(list.files(), read.csv)

Then carry one with code as above:

bind_rows(myReps, .id = "Rep") %>%
...etc
Question: Hi all,

I am performing a differential expression analysis with DESeq2, with these data:

- control, 2 replicates
- treatment, 2 replicates

So far there is still one obscure part of the manual for me: the *design* variable that you can set in many commands. I grasped the concept behind it but I am still struggling to understand how to use it properly. At the moment, I am using:

design = ~ condition

How would you set it, and why? Could someone write a couple of lines on how should I use that variable properly?

Any help appreciated!

Answer: In your case I would simply use the folowing, because the expression depends on only one factor which is the condition (either control or treatment).

condition = as.factor(c("control","treatment")
design = ~ condition

but if your replicates were not processed together, I would also take the batch effect into account.

batch = as.factor(c("rep1", "rep2")
design = ~ condition + batch

More examples on[ this tutorial][1].

[1]: http://www.bioconductor.org/help/workflows/rnaseqGene/#the-deseqdataset-sample-information-and-the-design-formula
Question: Hi All,

Again on the issue of marking duplicates...

If you have single end reads tools like Picard MarkDuplicates mark as duplicates reads sharing the same 5'end coordinates and ignore the 3' coordinate (see here). Is this the right thing to do?

I can see the point in looking at only 5' coordinates when the read length was quite a bit shorter then the fragment length (say 35-70 bp) and sequence quality dropped after ~50bp.

But now you can get reads up to 300bp, which is well into or above the range of the fragment length. So reads sharing the same 5'end but different 3'end might well be genuine independent fragments and MarkDuplicates would get them wrong. Am I missing something? Does anybody know of a tool that looks at both 5' and 3' ends for SE reads?

Answer: An update on this... I ended up writing a program to mark duplicates sharing start and end coordinates on the same strand and read group. The program is here MarkDupsByStartEnd.jar if anyone is interested.

@Kamil, I'm afraid samtools rmdup also look for 5'end only.
Question: [The following graph][1] shows an example cell trajectory analysis with monocle.

I have a rough understanding of how monocle constructs the pseudotemporal ordering of cells, using independent component analysis and minimum spanning trees (MST). However, there is one bit of the larger picture, which I'm still missing. How can monocle tell apart the most progenial cell (the source node; the beginning of the pseudotime) from the most differentiated cell (the sink node; the end of pseudotime)?

How does monocle know the direction of progress?

I can think of two possible type of answers to this question:

1. Pseudotime can in fact be reversed with respect to the real, biological time (i.e. beginning of pseudotime actually corresponds to the most differentiated cell, and the end of pseudotime corresponds to the most progenial cell).
2. There is something I don't understand about the biological significance of how monocle computes the minimum spanning tree, the associated diameter path, and the associated PQ tree -- and the process can, in fact, detect biologically correct root node, or it can detect it most of the time.

It is perhaps worth noting that another technique to establish cell trajectory, i.e. RNA velocity, has a very well understood method of establishing precedence, based on differential equations and mRNA biology.

[1]: https://www.nature.com/articles/nbt.2859/figures/2

Answer: Monocle does **not** know/determine the direction - it just finds the endpoints (and chooses one). As you can see from the [vignette][1] there is an argument called "root_state" which allows you to control where the pseudotime starts. There was however a recent paper which infered the direction - you can find it (and their method) [here][2].

[1]: http://cole-trapnell-lab.github.io/monocle-release/docs/#trajectory-step-3-order-cells-along-the-trajectory
[2]: https://www.nature.com/articles/s41586-018-0414-6
Question: how do I convert fa files to bed format?

Answer: You can't.

A Fa/Fasta file describes a sequence of DNA/Protein:

>Name
ATAGCTACGATTACGACGTACG
ATCGATCGATGCATCAGCTACT
AACTAGTCGATGATGCATACG...

A bed file describes some features mapped on a genome/sequence:

chr1 786 9879 gene1
chr2 486 979 gene2

The only thing you can do is saying that the BED contains only one feature: your sequence:

Name 0 1098 Name
Question: Hi everyone,

I used the factoextra package to make the biplot with top=20 contributed vars.
The code is:

**fviz_pca_biplot(pca_result1, title = "", col.var = "steelblue", select.var = list(contrib = 20))**

The result is very nice for me. However, I have a problem that I need the name of those 20 vars in a list format that I can use later.
Do you know any way to do that, except writing down the name of top=20 vars from the biplot graph?

Thank you,

Answer: Hi [@luongthang1908][1], following what is stated in the factoextra package documentation, you can get the top N contributing variables to both PC1 and PC2 by doing this:

pca_result1 = FactoMineR::PCA(X) # X is a dataframe/matrix with n rows (samples) and p cols (numeric values)
contrib_PC12 = pca_result1$var$contrib[,1:2] # get contributions of each variable to PC1 and PC2
eig_PC12 = pca_result1$eig[,1][1:2] # get the eigenvalues of both components
contrib_total = apply(contrib_PC12, 1, function(x) {sum(x*eig_PC12)/sum(eig_PC12)}) # calculate the total contribution
N = 20
names(sort(contrib_total, decreasing = TRUE)[1:N]) -> topN_PC12

Then `topN_PC12` contains the variables you are interested in.

[1]: https://www.biostars.org/u/56459/
Question: Hello everybody,

I am using the removebatcheffect function with simulated dataset using rnorm in order to test it.
So i have one experimental factor with two groups and one batch effect with 2 groups. The design is complet and balanced.
The correction seems to work well when i visualizing data on boxplot before and after correction but i dont understand why fold-changes between my two experimental groups before and after are the sames.

I tried the function using only the batch argument and also with the design argument but fold-changes are the sames each time.

Did anybody observe this using these function ?

Thank you.

Answer: I think i found a solution to my problem in looking removebatcheffect function in details . The explication come from the experimental design: while my batch effect is balanced through the experimental factor, so the batch estimated variance quantity that is substracted to the original data is balanced too through samples. Therefore means would not change even if values are not the same.
I thought i could have bigger foldchanges in removing the batch effect but that was wrong with this experimental design.

----------
Question: I have calculated XP-EHH and iHS scores for a set of snps using selscan. XP-EHH ranges from -0.75 to 0.9. What do extreme values show? In the original publication they plot log(P-value). I think that P-values was calculated to show differences between xp-ehh and iHS. Does anyone have experience with selection scores?

Answer: XP-EHH is a cross-population test for positive selective. It means that it detect SNPs that are under selection in one population but not in another. For example a SNP associated to resistance to malaria may be under selection in populations where malaria has been endemic but under neutral selection in other populations. It was developed because it is generally very difficult to identify signals of selection in a population, and comparing two populations may allow to identify weaker signals that are evident only after comparing with a closely related population.

The sign of the XP-EHH score indicates which of the two alleles is under selection, e.g. whether the ancestral or the derived allele. For practical reasons, people usually tend to ignore the sign and use the absolute XP-EHH score. This is because you may not always be sure about which SNP is ancestral in which population. Moreover taking the absolute score makes it easier to calculate mean by sliding windows.

In the publication they used -log(p-value), probably as a way to simplify the interpretation of the data in the plot. The value is usually generated by sorting the scores and taking the rank of them - e.g. see how I answered Zev in this discussion: https://www.biostars.org/p/94421/#163632 . The p-value is then converted to -log(p-value) to facilitate the interpretation. For example a p-value of 0.01 becomes -log(0.01) = 2, so you can say that all SNPs with a -log(p-value) higher than 2 are significantly selected in one population.
Question: I'm trying to learn the theory behind various steps in variant calling using GATK. Before alignment using BWA-MEM we first index the reference genome and this generates a set of files with the extensions

````chr13and17.fa.amb
chr13and17.fa.ann
chr13and17.fa.bwt
chr13and17.fa.pac
chr13and17.fa.sa````

where chr13and17.fa is the FASTA file containing the reference genome.

The next step in the pipeline is generating a .fai using samtools with the command:

````samtools faidx chr13and17.fa````

Followed by generating a .dict file using Picard:

`````java -jar picard.jar CreateSequenceDictionary
R=chr13and17.fa
O=chr13and17.dict`````

I want to know WHY we generate a .fai file and a .dict file despite also indexing the genome. In the samtools manual, the reason for creating a .fai file is specified as:

> Using an fai index file in conjunction with a FASTA/FASTQ file containing reference sequences enables efficient access to arbitrary regions within those reference sequences.

Isn't 'efficient access to arbitrary regions of the genome' also the aim of indexing? I understand the files themselves store different information in different, well, formats. But why all the different files though?

Answer: index for bwa-mem : burrow-wheeler transform index used to map the reads.

index fai : used by the tool to list the chromosome and quickly fetch a sequence from the fasta sequence

dict: list the chromosomes but also provides informations about the MD5 Sum of the fasta sequences (to be sure that you're using the same REF), the name of the organism(s), the names for aliases, the URL where we can retrieve the sequences, etc... this dict file will be inserted in/compared with the BAM and VCF headers
Question: Dear, I would like some help! I need to implement a computational infrastructure for my routine analysis of bioinformatics. I would like a suggestion on what would be a minimal or optimal configuration for tasks which can include: genome assembly and annotation, analysis rnaseq (transcriptome), then comparative genomics.

Thanks for attention.

EDIT: Post title changed to make it more informative by Ashutosh Pandey.

Answer: You should elaborate little more. See below:

1. Will it be just for your use? Or will it be part of a small research lab or some bioinformatics centre with multiple users?
2. Are you considering buying a multi CPU systems or cluster or you are talking about a desktop with good computing power and space (A few of those analyses can be done on desktop machines with descent RAM and memory)
3. Amount of storage capacities will depend on how much new data you are generating and analyzing. We have a md1000 systems that can take 8 X 500 Gb hard disks (4 TB). It is only used for the purpose of analysis and data storage for 6 months once it is out of the machine. We then need to back up those hard disks once every 6 months as our sequencing facility keeps generating new data.

Some related posts:

- https://www.biostars.org/p/11896/#11900
- https://www.biostars.org/p/64231/#64291
- https://www.biostars.org/p/63401/#63404
Question: Hi,

I have a problem with deleting columns in one of the networks in Cytoscape file.

Formerly the session file was big, with huge network from which I separated few small.
I also used Dynet app to compare the networks and this generated new columns.
Then I deleted majority of networks to create smaller file, which I assumed would be smaller on disk and will be loading faster (what wasn't the case but this is unresolved question).

So now I have six networks from which I want to delete some table columns to export whole file to html.
For one of the networks there are columns which are un-deletable. It is the case for a subset from other network in which I succeeded in deleting columns.
More curiously (and annoying) the unwanted coulmns (reamins from DYnet analysis) shows (and are un-deletable) after recreating the network from original one (which doesn't have this columns).

I thought that I could hide unwanted columns and in html-exported file they doesn't show but it does.

I'd be very grateful for help.

Answer: Ok, so I managed to do what I needed.

Maybe for someone it would be useful.

In short:

1. I unpacked the `cys` file and renamed to `zip`
2. In `tables` subdirectory I searched files with instances of "DAP" (the part of a columns' name, which I wanted to delete) and writeit to file `grep -rl 'DAP' > x `
3. I used `xargs` to delete files with 'DAP' hits `xargs -I{} rm -r {} < ~/x`

The file is still huge (and in the zip there are many not needed remains of deleted networks) but the troublesome colums disappeared.
Question: I'm working with pairwise whole genome alignment and retrieved some data from Ensembl, in the data there two types of symbols '-' and '.' . for example in this alignment:

http://www.ensembl.org/Homo_sapiens/Location/Compara_Alignments?align=714&db=core&g=ENSG00000184029&r=21%3A37951425-37973764

i think '-' represent gaps but i dont know what dots represent. in the example above there are dots in the vicinity of the sequence 'TTTAAAGAGAATATCCTTGTTCTTAGTAAATATGCCTT'

I appreciate if anybody could help me here

Answer: Answer by Ensembl Helpdesk:

> In this view:

> http://www.ensembl.org/Homo_sapiens/Location/Compara_Alignments?align=688--pan_troglodytes--17:62917953-63020080;db=core;r=17:64155265-64255264

> dashes ('-') represent gaps in the sequence of one genome relative to its counterpart in the pairwise view. You can think of them as indels (insertions/deletions when comparing the genomes).

> The dots are to represent that we were unable to find a match in that given region and this can be due to several region s such as the genome is fragmented (therefore regions are missing) or the match did not pass our threshold so we have not showed them in the alignment...

> This information should be in the help page of that view:

> http://www.ensembl.org/Help/View?id=202

> We will be updating it for the next release.
Question: I am working with certain codes in python and I wanted to know if a certain value of a variable can be accesssed outside the function without explicitly declaring it a global variable.
for example:

I was using this:

def func1():
v1 = "value1"
return v1
print ("This is the first function")

def func2():
v2 = "value2"
func1()
print ("This is the second function")

fo = input ("Enter 1/2\n")
if fo == '1':
func1()
elif fo == '2':
func2()
else:
print ("invalid")
print (v1)

It results me this error:

print (v1)
NameError: name 'v1' is not defined

Can I use this value of `'v1'` variable outside the function without defining it to b `'global'`??

Answer: You won't find many programming languages where all variables have global scope. If you want to access `v1` outside of `func1` then use `v = func()` in the main function and ensure `func1` returns `v1`.
Question: I posted this question on stackoverflow, but it did not get a response. It is a bioinformatics query, so perhaps this is a better forum: Is there are pythonic equivalent to the R ranges operation below?

In R (albeit longwinded):

Here is a test data.frame

df <- data.frame(
"CHR" = c(1,1,1,2,2),
"START" = c(100, 200, 300, 100, 400),
"STOP" = c(150,350,400,500,450)
)

First I make GRanges object:

gr <- GenomicRanges::GRanges(
seqnames = df$CHR,
ranges = IRanges(start = df$START, end = df$STOP)
)

Then I reduce the intervals to collapse into new granges object:

reduced <- reduce(gr)

Now append a new column to original dataframe which confirms which rows belong to the same contiguous 'chunk'.

subjectHits(findOverlaps(gr, reduced))

Output:

> df
CHR START STOP locus
1 1 100 150 1
2 1 200 350 2
3 1 300 400 2
4 2 100 500 3
5 2 400 450 3

How do I do this in Python?

Answer: https://github.com/biocore-ntnu/pyranges

import pyranges as pr
chromosomes = [1] * 3 + [2] * 2
starts = [100, 200, 300, 100, 400]
ends = [150, 350, 400, 500, 450]
gr = pr.PyRanges(chromosomes=chromosomes, starts=starts, ends=ends)
print(gr.cluster())
# +--------------+-----------+-----------+-----------+
# | Chromosome | Start | End | Cluster |
# | (int8) | (int32) | (int32) | (int64) |
# |--------------+-----------+-----------+-----------|
# | 1 | 100 | 150 | 1 |
# | 1 | 200 | 350 | 2 |
# | 1 | 300 | 400 | 2 |
# | 2 | 100 | 500 | 3 |
# | 2 | 400 | 450 | 3 |
# +--------------+-----------+-----------+-----------+

It will be out in 0.0.21. Thanks for the idea!
Question: I want to download the annotation file in gff3 format for the corresponding genome. As this fairly easy on the ncbi-webpage I don't find a possibility to do the same with efetch or the kind.

I hoped I could use something like this:

esearch -db nuccore -query "$genome_id" | efetch -format gff3 > "$path_data/$genome_id.gff"

Answer: I don't think NBCI offers GFF formatted files through efetch (yet).
Probably the best you can do is either indeed do it manually on the website (if you don't have many to do) or efetch genbank format and convert that to gff.

Otherwise, depending on the organism(s) you look for, there might be 'dedicated' databases that offer direct gff download.
Question: I used the SortSam option with GATK4 picard tools but it gives an error saying: htsjdk.samtools.util.RuntimeIOException: Problem writing temporary file file:///tmp/xxxxxx//sortingcollection.1885090825398510531.tmp. Try setting TMP_DIR to a file system with lots of space

I referred to this (https://www.biostars.org/p/42613/) post with the following command:

java -Xmx2g -Djava.io.tmpdir=`pwd`/tmp -jar SortSam.jar SORT_ORDER=coordinate INPUT=input.bam OUTPUT=output.sort TMP_DIR=`pwd`/tmp

However, I run GATK directly from the Linux terminal with the command:

./gatk SortSam
Also, there is no SortSam.jar file in the directory. Only a gatk-package-4.0.11.0-local.jar.

How do I modify the syntax of the above command to set the temporary directory ?

Answer: I got it using

./gatk --java-options "-Djava.io.tmpdir=/custom file path" SortSam -I=samplealn.sam -O=sortedsamplealn.bam -SO=coordinate --TMP_DIR=/custom file path
Question: What does this mean: non-redundant template structures

Answer: template structures: known [protein|rna|other biomolecule] structures that will be used in homology-modeling and threading/fold recognition to generate a structure for your input sequence based on sequence similarity

non-redundant template structures: a set of all known structures that are deemed to be unique in the set of all known template structures. Uniqueness is determined according to some similarity measure (e.g. backbone RMSD of superposed 3D structures, sequence similarity of primary structures, topology differences of secondary structures etc). Non-redundancy allows you to limit your search space, otherwise you'd try using the same (or very similar) template structure all over again, which increases overall runtime and complexity.
Question: Does anyone know where I can find a table that gives the dbSNP build numbers alongside their release date?

I know dates are available for each release in the summary information (e.g. http://www.ncbi.nlm.nih.gov/SNP/snp_summary.cgi?view+summary=view+summary&build_id=129), but I'm wondering if there's a full list somewhere.

Thanks in advance.

Answer: http://www.ncbi.nlm.nih.gov/projects/SNP/buildhistory.cgi