cannin/biostars_qa huggingface dataset:
Question: How can I merge 2 vcf files, each with multiple samples, the so the final file has all the samples but the intersection of the SNPs. How about union on the SNPs?

Answer: There are multiple questions in your post but I'll make some answers that I believe will answer all of them.

Firstly, I recommend that you normalise your VCFs so that everything else that you do will work as planned. 'Normalising' a VCF usally means splitting multi-allelic calls, left-aligning indels, and ensuring that the alleles specified as REF in your VCF are also present in the reference genome. You can do this with (for each VCF file):

#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
bcftools norm -m-any MyVariants1.VCF | bcftools norm -Ob --check-ref w -f /ReferenceMaterial/1000Genomes/human_g1k_v37.fasta > MyVariants1.norm.bcf ;
bcftools norm -m-any MyVariants2.VCF | bcftools norm -Ob --check-ref w -f /ReferenceMaterial/1000Genomes/human_g1k_v37.fasta > MyVariants2.norm.bcf ;

bcftools index MyVariants1.norm.bcf ;
bcftools index MyVariants2.norm.bcf ;

You can then merge these together with:

bcftools merge -Ob -m none MyVariants1.norm.bcf MyVariants2.norm.bcf > MyVariants.norm.merge.bcf ;

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

For the intersection, etc, you may additionally consider setting a unique value for each variant in the ID column. The following line will eliminate whatever is already set for the ID field and then set a new value:

bcftools annotate -Ob -x 'ID' -I +'%CHROM:%POS:%REF:%ALT'

You can then do various operations using `bcftools isec`. There are examples given in the program when you just type `bcftools isec`:

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

Hope that this helps.

I have done intersections and other complex joins before, but by using AWK indexed arrays (can use Python, too) after having set the ID field to a unique value, like I did in the code above. Sometimes, *BCFtools*, *VCFtools*, etc., do not behave as you expect.

Kevin
Question: Hi everyone,

Wondering if someone can provide me with some guidance.

I have previously sequenced 4 skin cancers using 10X chemistries and I would like to combine them into one dataset. My research question is to look at cancer stem cell populations, so I will need sensitivity.

I have done a load of reading and I am still confused as to what the correct method is.

CellRanger Aggregate causes a massive downsampling in 2 out of 4 samples. I believe this is becuase the number of cells in each sample varies (~600 - ~2500 cells). If I proceed with the aggregated dataset, I will be losing too much sequencing depth which will hinder my search for the cancer stem cells. Another option is to use CellRanger Aggr's 'normalize=none' argument, which does not normalize the average read depth per cell between samples before merging. From here, could I then go onto use scTransform or Seurat's 'NormalizeData' arguments?
I have performed QC for each sample dataset and they all need different filtering parameters (nFeature_RNA), so the only problem with this method would be giving the combined dataset one filter, which could potentially include bad data for some of the samples.

Another option that I could do is filter and normalize each sample individually and then use Seurat's 'merge' function. However, by default, this uses the raw count matrices which will erase the normalized data matrices. There is an option to merge the normalized data matrices as well as the raw data matrices (merge.data = TRUE) but this should only be done if the same normalization approach was applied to all the objects. Does this mean the same normalization method such as 'LogNormalize' or 'CLR', or does it mean the same filters?

I am sorry if this is a disgruntled 'question', I am just looking for the best method to combine my 4 datasets without loosing too much data and then continue with downstream analyses.

Thank you so much in advance!!

Alex

Answer: Avoid CellRanger Aggregate.

Using Seurat `merge` with `merge.data = TRUE` is fine. It just means the same normalization method, not the same filters. You *should* define your filters for each sample separately.
Question: Hi,

I have a list of transcripts, such as: ENST00000416191. I want to find the Ensembl IDs for the genes they are associated with. Would I get identical results, whether I query GRCh38 or GRCh37? Note that I am not interested in their coordinates, only the IDs.

Thank you

Answer: Hi jiab,
to give you an idea of the differences here is a little bash script, which gives you a list of ENSTs mapping to different ENSGs in grch37 and grch38:

#!/bin/bash
#download of all ENST mapped to ENSG from Ensembls grch37 biomart
wget -O mapping_37.txt 'http://grch37.ensembl.org/biomart/martservice?query='

#download of all ENST mapped to ENSG from Ensembls grch38 biomart
wget -O mapping_38.txt 'http://www.ensembl.org/biomart/martservice?query= '

#join the two mappings on ENST
join <(sort mapping_37.txt) <(sort mapping_38.txt) >joined_on_ENST.txt

#get just the lines, where ENSGs differ for grch37 and grch38 for the same ENST
awk '$2!=$3 {print $1, $2, $3}' joined_on_ENST.txt >different_ENSG_37_38.txt

If you have a look at the number of lines in each file:

wc -l *
1421 different_ENSG_37_38.txt
190566 joined_on_ENST.txt
215170 mapping_37.txt
218207 mapping_38.txt
It shows you that the mappings in 37 and 38 share about 190k ENST having 25k (37) and 28k (38) unique ENSTs and 1421 ENSTs mapping to different ENSG in both builds.
Out of about 243k different ENSTs (190566 +25k+28k) you could search your ENSG for, in the different versions, you would not get the same result for both builds (no ENSG in one build, or different ENSG) for about 54k queries (25k+28k+1421).
Question: Hello,

I am trying to use the pysam package in python to do some processing of a large BAM file, which is a BWA output.
The header of my file looks like that:

$ samtools view -H myFile.bam
@HD VN:1.3 SO:coordinate
@SQ SN:1 LN:...
@SQ SN:2 LN:...
@SQ SN:3 LN:...
@SQ SN:4 LN:...
@SQ SN:5 LN:...
@SQ SN:6 LN:...
@SQ SN:7 LN:...
@SQ SN:8 LN:...
@SQ SN:9 LN:...
@SQ SN:10 LN:...
@SQ SN:10000001 LN:...
@PG ...

and I have the following alignment record within this file (just an example):

REC-XYZ-1 147 10000001 1 36 9S91M = 3601 3511 NAAAAACTCTAGTTGCTCGTGGGTGTCAACTTTTGATTATATGCTGTTATGTTTGTGCAAACTCAACTGAAATGTGAACTCAAGTGGTGATGCTATGAAC DBCC@:CDCCDDDDDDFFHGHIIIIEGHEHIIGIJIJIEHCGGFIDHHCHJJJJIIIIHHDIFGIFJIIIGIFJJJGHGEIIGIGGIHDFHHFFFDF@@@ NM:i:0 MD:Z:91 AS:i:91 XS:i:99
(notice the '10000001' in the 3rd field).

Using pysam, I do the following:

bam_file = pysam.AlignmentFile("myFile.bam", "rb")
it = bam_file.fetch("10000001")
for x in range(20):
aln = next(it)
print(aln.reference_id)

10
10
10
.
.
.

surprisingly, although the real reference id is '10000001', pysam sees it as '10'. For example, the same record I showed above looks like this:

REC-XYZ-1	147	10	0	36	9S91M	10	3600	91	NAAAAACTCTAGTTGCTCGTGGGTGTCAACTTTTGATTATATGCTGTTATGTTTGTGCAAACTCAACTGAAATGTGAACTCAAGTGGTGATGCTATGAAC	array('B', [35, 33, 34, 34, 31, 25, 34, 35, 34, 34, 35, 35, 35, 35, 35, 35, 37, 37, 39, 38, 39, 40, 40, 40, 40, 36, 38, 39, 36, 39, 40, 40, 38, 40, 41, 40, 41, 40, 36, 39, 34, 38, 38, 37, 40, 35, 39, 39, 34, 39, 41, 41, 41, 41, 40, 40, 40, 40, 39, 39, 35, 40, 37, 38, 40, 37, 41, 40, 40, 40, 38, 40, 37, 41, 41, 41, 38, 39, 38, 36, 40, 40, 38, 40, 38, 38, 40, 39, 35, 37, 39, 39, 37, 37, 37, 35, 37, 31, 31, 31])	[('NM', 0), ('MD', '91'), ('AS', 91), ('XS', 99)]
(notice the '10' in the 3rd field)

The same thing happens for the RNEXT field (which I get through aln.next_reference_id), and this is really bad because I cannot differentiate between chromosome 10 and chromosome 10000001.
Does anyone know why this is happening? Is this a bug in pysam? I would be quite amazed to find out that it can only take two-digit values as reference id (nothing about that in the [documentation][1]).
Any ideas?

\* I am running on linux, with python 2.7 and the latest version of pysam, using IPython notebook.
[1]: http://pysam.readthedocs.org/en/latest/api.html

Answer: If you look at the [read the docs](http://pysam.readthedocs.org/en/latest/api.html#sam-bam-files) page for pysam, you'll see the following warning regarding `reference_id`:

> This field contains the index of the reference sequence in the sequence dictionary. To obtain the name of the reference sequence, use `pysam.AlignmentFile.getrname()`

In other words, pysam isn't confusing anything, it just happens that many of the IDs happen to correspond to chromosome names in your case (though they'll be the **wrong** chromosome names).
Question: I'm new to the field of metagenomics and I have been reading research work done on binning in metagenomics.

I have been reading this [paper][1] where the authors discuss about different metagenome binning approaches. The last approach is **Differential abundance (DA)**. They explain the starting point of this method as *differential coverage patterns across multiple samples where population changed in abundance over time*. Can someone kindly explain to me what is meant by **differential abundance**? Is it the abundance of genomes in different samples?

Sorry if this is a naive question. I would really appreciate an explanation. :) Thank you!

[1]: https://microbiomejournal.biomedcentral.com/articles/10.1186/s40168-016-0154-5

Answer: I believe what they're referring to is the relative change in genome abundances between samples where the genome distributions change relative to time. The pattern they're referring to would be the common distribution shift between samples, which would allow clustering of samples based on their common genome enrichment.

A toy example of this would be if there were 4 genomes sequenced across 100 samples and were uniform across all 4 at the start of the time series. At the end of the time series, if there were changes in abundances, it would be possible to use DA to group the samples that had the same enrichment in their genome fractions. If 25 of the samples all enriched for genome 1, they could be considered more similar than the remaining 75 which down selected genome 1.

I think in this reference though they are talking about contigs and not full genomes, but the idea is the same.
Question: I am trying to calculate the number of paralogs for a few genes in different species in Ensembl. I am wondering if there is any tool which can do it automatically. Thank you.

Answer: As from your question it looks you like to use Python for this, you could also opt to use the Ensembl REST API.

It's easy to get all orthologues for a gene, for example for the human ABCD1 gene:

http://beta.rest.ensembl.org/homology/symbol/human/ABCD1?content-type=application/json;format=condensed;type=orthologues

or all paralogues for the same gene:

http://beta.rest.ensembl.org/homology/symbol/human/ABCD1?content-type=application/json;format=condensed;type=paralogues

So, using these REST statements, I think it should be quite easy for you to start out with a particular gene in e.g. human, find the orthologues in the other Ensembl species and then get the number of paralogues for those genes.

How you can use REST statements in Python code you can find in the REST documentation for the statement in question.

Hope this helps.
Question: I've noticed that one of my .sorted.bam files contains a read that does not show up on SAMtools tview. The bit flag associated with that read suggests it's an optical or PCR duplicate. Does anybody know if SAMtools tview hides these kind of reads by default or is something else potentially happening?

Answer: In the old versions of samtools when pileup still existed, the default settings for pileup was to filter reads with bitwise flag 0X704.

So for pileup generation the following reads were not considered at all from the bam files:

1) 0x0400 (aka 1024 or "d") duplicate

2) 0x0200 (aka 512 or "f") failed QC

3) 0x0100 (aka 256 or "s") non primary alignment

4) 0x0004 (aka 4 or "u") unmapped

I assume tview may be following the same filtering strategy and not showing reads that are PCR duplicates, unmapped and non-primary alignments. This makes sense too because sometimes PCR duplicates may create visualization problem if there are large number of duplicates reads at a given position.
Question: I am trying to use the [UniProtJAPI for Java][1] to get the protein sequence based on the entry name and later identifying a consensus sequence. I'm working on NetBeans 8.2 and UniProt JAPI 1.0.14

I have a list of entry names along with the peptides found by mass spectrometry, a short example (not actual data) is:

> MAOM_YEAST_R.LATYGGD.K

**MAOM_YEAST** is my key to the full sequence, after separating it from the partial sequence I want to use this to get the UniProt ID and from there the corresponding full protein sequence. This full protein sequence is to be used to extend the partial sequence to find possible matches outside of the range of the partial sequence. So to quickly sum up the steps:

> MOAM_YEAST > find UniProtID (P36013) > find sequence > locate partial
> sequence > extend the partial sequence with ~5 amino acids in both directions > search
> for a consensus sequence

I have consulted the **UniProt JAPI documentation** (included in the download), but especially **uk.ac.ebi.kraken.interfaces.uniprot / UniProtId** is where the confusion starts. To cite from the documentation:

> How to work with this Interface
>
> The standard way of retrieving this data type
>
> The standard way of setting this data type
>
> UniProtEntry entry = getEntryFromParserOrAPI();
>
> entry.setUniProtId(DefaultUniProtFactory.getInstance().buildUniProtId("CYC_HUMAN"));
>
>
> UniProtId id = entry.getUniProtId();

However, I get errors for the UniProtEntry entry not being loaded properly due to the **getEntryFromParserOrAPI()** not working (**don't have the exact error at the moment, will post it ASAP**). This method seems like the ideal way to perform the action I want, replacing "CYC_HUMAN" with another name to get the proper entry. If I understand the documentation correctly a UniProtEntry should be able to get the UniProt ID based on an entry name like "CYC_HUMAN" using getUniProtID() and getSequence() could be used for the sequence.

My questions are:

**1) Does anyone know how to go from an entry name like MOAM_YEAST or CYC_HUMAN to the corresponding UniProt ID and perhaps from there to the sequence?**

**2) Does anyone have a solution for the suggested code from the documentation to get it working?**

Much thanks

[1]: https://www.ebi.ac.uk/uniprot/japi/

Answer: 1) Does anyone know how to go from an entry name like MOAM_YEAST or CYC_HUMAN to the corresponding UniProt ID and perhaps from there to the sequence?

$ curl -sL "http://www.uniprot.org/uniprot/MAOM_YEAST.xml" | \
xmllint --xpath '//*[local-name()="sequence"]/text()' - | \
tr -d '\n' |\
awk '{S="LATYGGD";x=5;i=index($0,S);print substr($0,i-x,length(S)+2*x);}

SIECRLATYGGDKDVDY
Question: Hello,

I am trying to use BUSCO to assess the quality of a transcriptome assembled with Trinity. My run fails after ~2h running, when the HMMER phase begins. I couldn't find any support options for BUSCO.

Is there any way to get support for BUSCO? I couldn't find anything on line or on the documentation.

This is my command and programs versions:

```
module load Python/3.2.1
module load BLAST+/2.2.31
module load HMMER/3.1b2
module load EMBOSS/6.5.7

python3 ~/MyPrograms/BUSCO_v1.1b1/BUSCO_v1.1b1.py -in ../Trinity-GG.fasta -o GG \
-l ~/MyPrograms/BUSCO_v1.1b1/vertebrata -m trans -c 4 -f
```

After ~2 hours running, this is the error I get:

```
*** Running HMMER to confirm transcript orthology ***
Traceback (most recent call last):
File "/.../MyPrograms/BUSCO_v1.1b1/BUSCO_v1.1b1.py", line 594, in
name=i[:-7];group=transdic[name]
KeyError: 'TRINITY_'
```

Any help would be appreciated

Thanks

Answer: Hi

Thank you both for your suggestions. I was able to solve this without using docker.

I tried version 1.1 but I got the same error.

The newer versions of Trinity dropped the troublesome pipe (|) from their headers, so most compatibility issues are now solved. However, the problem was indeed in the header: it won't accept Trinity's header for genome-guided assemblies:

TRINITY_GG_XXXXX

Deleting _GG worked for me:

sed s/TRINITY_GG_/TRINITY_/ Trinity-GG.fasta > Trinity-GGbusco.fasta

I had no problems with running BUSCO with DeNovo assemblies, which have the header TRINITY_DNXXXXXXXXX

Thanks
Question: Hi everyone,

I am mid-way through my university project, and I have run into an issue. I have a long list of around 1000 proteins that I wanted to analyse in STRING, however, my list is too large. I decided to try and utilise Cytoscape (and downloaded the stringApp), but the networks generated are still very messy. I've attached a screenshot: https://i.stack.imgur.com/KUIE4.jpg. Is there any way to improve the presentation of the network by downloading any Cytoscape apps or by tweaking the settings? I am very much a novice when it comes to Cytoscape, even after watching quite a few tutorials online.

Thanks in advance

Answer: It would be very difficulty to see individual proteins and their interactions when the number is 1000. In such cases the best you can do is play with different layouts in Cytoscape to get a meaningful global representation. I have made networks with 15000+ proteins in Cytoscape, but those are not meant for looking into individual proteins and their interactions.

Here is an example of all `c.1.` folds in SCOP and their relationships, which is comparable in number to what you have.

![enter image description here][1]

As you were already told above, adjust the stringency of your criteria to get most reliable interactions for your protein(s) of interest, and that should reduce the number of proteins and be easier to visualize. Either way you will have to play with layouts, or else any network will end up looking like a blob in your image.

[1]: https://i.ibb.co/vqnRV1C/SCOP-175-number34-Prob-e06-P10-c-1-fold.png
Question: Hello. I have a dataset of cases and a dataset of controls which have been genotyped on different platforms. So, for a particular SNP, it can happen that I have the same allele coding both in cases and in controls; or it can happen that I have in controls the reverse of cases (e.g.: A C in cases and T G or G T in controls).
I want to merge the files of cases and controls, but the allele coding should be consistent: which procedure can I follow? Please, let me know. Thanks!

Answer: One possible workflow is described at http://www.cog-genomics.org/plink/1.9/data#merge3 .
Question: Hi,

I'd like to automatically download the tar files from the following links, does anyone knows how to do that?

http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39977

http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36375

Kind of like using wget...

Thanks!

Answer: Just a note that the GEOquery bioconductor package does this for a living, if you happen to be targeting R for further processing/analysis. See the help for getGEOSuppFiles() and getGEO().
Question: Dear Friends,
I want to download the entire human gene list with the information about their chromosomal location, i.e. that specific gene is present in which chromosome and which arm (p or q). Is there any specific databases to give such an information, please guide me.
Thanking in advance for the replies.

Answer: Just use Biomart. Using the Ensembl annotation, you can get a list of Ensembl IDs, the associated gene names, chromosome, position, and band (such as p22.1). You can also filter the results for the exact type of gene (known, novel, etc.) you want.

Edit: FYI, if you just want the actual coordinates instead, just download a GTF file or follow one of the other answers.
Question: There are different terms for the quality/status of a genome assembly. Like "draft genome", "finished genome" or "high quality genome".

For the case of bacteria I am interested in any standard/"official" definition of quality standards. I basically understand the meaning of the definitions above. But I would like to have a source where I can find an overview of this kind of definitions (what else is commonly used?) and what they are supposed to describe. As already stated, I am mainly dealing with bacteria, so the definitions would ideally apply to bacteria.

Answer: Check the Sequence ontology. They have a controlled vocabulary for the qualities of the genomes.

For the term '[whole_genome_sequence_status][1]' you find definitions for the following children:

```
noncontiguous_finished (SO:0001490)
high_quality_draft (SO:0001487)
standard_draft (SO:0001486)
improved_high_quality_draft (SO:0001488)
finished_genome (SO:0001491)
annotation_directed_improved_draft (SO:0001489)
```

[1]: http://www.sequenceontology.org/browser/current_svn/term/SO:0001499
Question: Hi guys,
is it possible that there exist tumor patients without mutations in driver cancer genes?

Kind regards,

B

Answer: Obviously, certain viruses can initiate oncogenesis by introducing their own drivers into the host cell (rather than mutating host-cell drivers).
But also, it's very likely that solid tumours may lack mutations in known driver genes, since these gene lists are continually being redefined as more studies are performed.

... and look into medulloblastoma

... and not all mutations in driver genes are driver mutations
Question: Hi,

I have sequence from grapevines and I'm trying to see what viruses there are. I was able to assemble most of a viral genome (15000 out of 18000 bp in one contig). I'm trying to estimate the coverage. Here's a breakdown of what I did:

raw reads -> trim adapters -> map to grape reference and remove mapped reads -> assemble trimmed, unmapped reads.

To estimate the coverage, I used the virus that had the largest contig (grapevine leafroll associated virus 3) and I mapped the trimmed unmapped reads to it. I started with about 7 million reads and 1.3 million of them mapped to the genome. The average read length was 50 bp and the total genome size was 18Kbp. Using the equation presented in this: http://res.illumina.com/documents/products/technotes/technote_coverage_calculation.pdf I get

Coverage = Length of read * number of reads / haploid genome length

Coverage = 50 * 1.3x10^6 / 1.8x10^4

Coverage = 3611x? Could that be right?

Answer: For viruses, we often get that kind of coverage.

To sanity check it, I would open up the alignment in a viewer like IGV and see if it looks about right.

RNA viruses have really uneven coverage so you should expect to have regions of really high and (alas) no coverage. I don't know about DNA viruses. It may be the RNA secondary structure that makes it uneven.
Question: Hi guys,

I am not sure, in terms of sequencing what is referred to as color space and nucleotide space. What is a difference between a normal reference genome and a color space genome. Does it has to do something with base calling, that the original data is retained without the conversions to textual nucleotide data!!!! What we use in normal chip-seq and rna-seq analysis.

There's a kind of similar question here but the answer is missing the link.

Thanks

Answer: Here is a tutorial I wrote on this topic:

- https://www.biostars.org/p/43855/