cannin/biostars_qa huggingface dataset:
Question: Hi all,
I am following the "**best practice**" suggested by broad institute to call variants from whole exome sequencing data. Currently, I am using **Mutect2** to call variants from **tumor** sample and **normal** sample based on latest reference genome **GRCh38.**
But, I don't have *interval list* to use **-L** option.
Is there any way to generate interval list from exome sample which I have? or Is there any default interval list for exome data?
Thank You

Answer: > Is there any default interval list for exome data?

something like:

$ curl -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg38/database/refGene.txt.gz" |\
gunzip -c | cut -f 3,5,6 | sort -t $'\t' -k1,1 -k2,2n | bedtools merge -i - > exome.bed
Question: Hi,

You please imagine I have 23 differenrially expressed genes between tumour (26 samples) and normal (30 samples) tissues.

I have calculated the Pearson correlation values between the expression of these 23 differentially expressed genes in individual samples (56 samples) and the mean of expression of this gene in tumour `(cor1)` and normal `(cor2)` like below

> head(pred)
cor1 cor2 scores samples
cl1 0.8824150 0.7623499 0.06003256 tomour
cl1.1 0.7928945 0.6831890 0.05485271 tomour
cl1.2 0.8434347 0.7933102 0.02506226 tomour
cl1.3 0.7174795 0.5120342 0.10272268 tomour
cl1.4 0.6571185 0.6084326 0.02434294 tomour
cl1.5 0.5273702 0.6716480 -0.07213889 tomour

> > dim(pred)
[1] 56 4
>

And I have calculated a score like so

(cor1 – cor2)/2 (range: −1 to +1)

+1 says this is tomour sample and -1 says this is a normail sample

With such tables how I can represent a correlation plot like this picture? in this plot nicely samples were divided by scores and correlation

![enter image description here][1]

By this code

b <- ggplot(pred1, aes(x = tomour, y =normal))

b + geom_point(aes(color =samples, shape = samples))+
geom_smooth(aes(color = samples, fill = samples), method = "lm") +
geom_rug(aes(color =samples)) +
scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))+
scale_fill_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))

I have plotted

![enter image description here][2]

But this is not similar to what I want at all

[1]: https://i.ibb.co/W5Gkdtn/Untitled.jpg
[2]: https://i.ibb.co/fCZs5j0/Rplot03.jpg

Answer: use geom_jitter from ggplot

library(ggplot2)
library(cowplot)

dat <- data.frame(
sample=1:10,
score=c(
rnorm(n=40,mean = 0.7,sd=0.1),
rnorm(n=40,mean=0,sd=0.1),
rnorm(n=20,mean=-0.5,sd=0.1)),
diagnoses=c(
rep("ADC",40),
rep("SCC",40),
rep("Other NSCLC",20))
)

head(dat)
# sample score diagnoses
# 1 1 0.5154322 ADC
# 2 2 0.7057089 ADC
# 3 3 0.6860727 ADC
# 4 4 0.5880083 ADC
# 5 5 0.7088662 ADC
# 6 6 0.6122866 ADC

# plot
ggplot(dat,aes(x=1,y=score,color=diagnoses))+
geom_jitter(size=2) +
geom_hline(yintercept = 0) +
geom_hline(yintercept = 0.2,linetype=2) +
geom_hline(yintercept = -0.2,linetype=2) +
theme(axis.title.x=element_blank(), # remove x-axis label
axis.text.x=element_blank(),
axis.ticks.x=element_blank())

![enter image description here][1]

[1]: https://i.imgur.com/6fdsWXP.png
Question: Hi All,

I would like to ask that is there an accurate tool/script for deinterleaving fastq files? I managed a 'bam to fastq' conversion according to http://gatkforums.broadinstitute.org/discussion/2908/howto-revert-a-bam-file-to-fastq-format and I would like to deinterleave the output fastq.

Many thanks!

Answer: In terms of performance, I have yet to see any program written in Perl, Python, C++, or C outperform [this combination of shell tools][1]. Other tools may be more convenient--and often that is more than a good enough reason to use them--but if you need blistering speed then there is no substitute for native built-in commands.

To summarize, the solution (using paste and process substitution) looks something like this.

```
paste - - - - - - - - < reads-int.fq \
| tee >(cut -f 1-4 | tr "\t" "\n" > reads-1.fq) \
| cut -f 5-8 | tr "\t" "\n" > reads-2.fq
```

[1]: https://gist.github.com/nathanhaigh/3521724
Question: Hi all,

I have three VCFs, a child (male), a father and a mother, and I would like to extract de novo variants in the child.

All three samples were called separately, however, using the same GATK pipeline. I ran rtg tool to try to extract the de novo variants following this command:

rtg mendelian -t /path/to/referencegenome.sdf --input trio.vcf.gz --lenient --output-inconsistent trio-non-mendelian.vcf.gz

According to the rtg manual:

>Records where the presence of missing
values makes the Mendelian consistency undecidable contain MCU INFO annotations in the annotated output VCF.
The following examples illustrate some consistent, undecidable, and inconsistent calls in the presence of missing
values:

Below are the headers related to rtg in my output vcf:

##INFO=

##INFO=

Below are some annotated info in my VCF. I'm having a difficulty in understanding the output and would hightly appreciate if someone could explain how to interpret the resulting annotations.

INFO	FORMAT	Child	Father	Mother
;MCV=Child:0/0+1/1->0/0 GT:AD:DP:GQ:PGT:PID:PL . . 1/1:0,2:2:6:1|1:10433_A_ACCCTTAACCCCTAAC:90,6,0
;MCV=Child:1/1+0/0->1/1 GT:AD:DP:GQ:PL 1/1:0,5:5:15:157,15,0 1/1:0,3:3:9:93,9,0 .
;MCV=Child:0/0+0/0->0/1 GT:AD:DP:GQ:PL 0/1:4,5:9:72:99,0,72 . .

Answer: Hello,

due to the use of lenient flag, you see that missing spots are assumed as 0/0 (which is mostly correct). In case you have a gVCF file, you can try repeating the analysis on it, and then all three samples should have genotypes assigned, and then the INFO annotation will fully make sense.
Question: Hello,
I am having a hard time opening a very large bedgraph. I have been suggested to split my bam file by chromosome with chrom-bed.py but it didn't work.
Is there any other alternative?
Thanks,
GP.

Answer: samtools view in.bam chr1 -b > out.bam

Use `-b` to output bam format
Question: Dear all,

I have FASTQ files and on start of my read I have 7 nucleotides tag - I would like to extract reads with this specific tag.

I would like to search in first 15 nucleotides of my reads, if match - extract this read to new fastq files.

Thank you for any ideas or help.

Answer: You can also do this with BBDuk:

bbduk.sh -Xmx1g in=reads.fq out=filtered.fq k=7 mm=f rcomp=f restrictleft=15 literal=ACGTACG

If you want multiple tags, you can list them separated by commas; and if you want to allow 1bp mismatch, set hdist=1.
Question: I found that GEO holds "processed sequence data files", while SRA holds "raw sequence data files". In which way processed?

I am interested in rna-seq bacterial data, is GEO right resource for me?

Answer: Depends on your needs. GEO contains data relevant to a particular experiment and in most cases I believe represent processed data (e.g. have been run through one or more analysis steps, such as trimming, alignment to a reference, R/BioC, etc). The example RNA-Seq from GEO is illustrative:

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

If you check under the BioSamples you'll find a Cufflinks file, a SAM file, BEDGRAPH, etc. These would all be bound to a specific assembly version and whatever comes along with that (annotation, etc).

SRA however contains raw sequence data for an experiment, if you want to download and re-analyze the data on your own. They are generally all tied together via a common BioProject ID.

**EDIT**: In that last sentence, by 'they' I mean any data relevant to the BioProject (BioSamples, SRA, assemblies, etc).
Question: I'm looking for a free software tool, R-package or Matlab script to create interactive 3D scatter plots of biological data (gene and protein expression data, or metabolite abundances), which allows the users to add biological annotations to the data points (representing biomolecules) and display them when a user hovers the mouse over them. I found many scatter plot tools so far, but adding interactive annotations does not seem to be supported in most cases.

Many thanks for your help.

Answer: You can try the R-package [vrmlgen][1] which provides annotated 3D scatter plots in VRML and LiveGraphics3D format (see also page 7 of the [documentation][2]). Works also with large numbers of data points.

[1]: http://cran.r-project.org/web/packages/vrmlgen/index.html
[2]: http://www.jstatsoft.org/v36/i08/paper
Question: Hello, I am a trying to filter some vcf files based off a set of rs#s (`$FILE` corresponds to chromosome #, so for now that is just 21, 22, and Y)

vcftools --gzvcf "sorted_AltaiNea.hg19_1000g."$FILE".mod.vcf.gz" --snps 330k.txt --out "filtered_AltaiNea.hg19_1000g."$FILE"_"

I run that code on the cluster, but then all I get out are

filtered_AltaiNea.hg19_1000g.21_.log
filtered_AltaiNea.hg19_1000g.22_.log
filtered_AltaiNea.hg19_1000g.Y_.log

Which contain text as such:

VCFtools - v0.1.11
(C) Adam Auton 2009

Parameters as interpreted:
--gzvcf sorted_AltaiNea.hg19_1000g.21.mod.vcf.gz
--out filtered_AltaiNea.hg19_1000g.21_
--snps 330k.txt

Using zlib version: 1.2.3
Versions of zlib >= 1.2.4 will be *much* faster when reading zipped VCF files.
Index file is older than variant file. Will regenerate.
Building new index file.
Scanning Chromosome: 21
Warning - file contains entries with the same position. These entries will be processed separately.

Writing Index file.
File contains 35104060 entries and 1 individuals.
Applying Required Filters.
Keeping sites by user-supplied list
After filtering, kept 1 out of 1 Individuals
After filtering, kept 5364 out of a possible 35104060 Sites
Run Time = 176.00 seconds

I do not have any vcf output, just those log files. Is there something I am doing wrong? Thanks!

Answer: You should add the `--recode` flag to produce a new vcf file.

See example 2 on http://vcftools.sourceforge.net/man_latest.html
Question: If you have a closed bacterial genome or plasmid, what are the general rules where to put the start. With start I mean where the base "1" of the finished sequence is (or you could also say with which gene the sequence starts).

I heard some rules already like putting the OriC at the beginning.

Does anyone have a good paper/source for that?

Thanks!

Answer: There is no general rule, but please avoid to linearize your circular sequence within a coding sequence. There are many plasmids in Genbank where interesting genes like antibiotic resistance determinates have been cut into two halves. Its always a pain to work with such sequences.

In most Eubacterial chromosomes the origin of replication is located between the highly conserved genes rpmH and dnaA. For Firmicutes it is a good choice to place the start codon of dnaA on position 500. In some species there may exist isolates where the orientation of the oriC region is reversed with respect to the majority of the chromosome. With such an isolate it may be appropriate to make rpmH the first gene on the linearized chromosome and dnaA the last.
Question: Hello all!

I have used OrthoFinder to identify orthogroups among five species. Part of the output is a rooted phylogenetic tree, but I think it is not ultrametric?

I would like to use the tree as part of the input for a CAFE analysis. The CAFE tutorial described using the program r8s to make a tree ultrametric, but I don't know the number of sites in the alignment used to estimate the species tree. TimeTree has approximate divergence time for some of the species in my tree, which I can use to calibrate the tree?

Any suggestions on how I can create an ultrametric phylogenetic tree would be greatly appreciated!

Nancy

Answer: If you are familiar with Python you can convert a rooted tree into ultrametric topology (all leaves will have the same distance to root) using [ETE Toolkit 2-3][1].

from ete2 import Tree

t = Tree('(A:1,(B:1,(C:1,D:1):0.5):0.5);')
# Root tree at B node.
t.set_outgroup('B')

print('Non-ultrametric rooted tree:')
print(t.write())

t.convert_to_ultrametric()
print('Ultrametric rooted tree:')
print(t.write())

Output:

Non-ultrametric rooted tree:
(B:0.5,((C:1,D:1)1:0.5,A:1.5)1:0.5);

Ultrametric rooted tree:
(B:2,((C:0.666667,D:0.666667)1:0.666667,A:1.33333)1:0.666667);

[1]: http://etetoolkit.org
Question: Hi,

I want to print a tree with selected attributes using ete3 in python using a Newick (with NHX annotations) file e.g. :

(dog[&&NHX:Genes=3:Identical=3:Duplicated=0:Lost=0:Novel=0:Singleton=0], ((mouse2[&&NHX:Genes=19:Identical=19:Duplicated=0:Lost=0:Novel=0:Singleton=0], mouse[&&NHX:Genes=21:Identical=17:Duplicated=4:Lost=0:Novel=0:Singleton=0])mouse_mouse2[&&NHX:Genes=19:Identical=9:Duplicated=0:Lost=1:Novel=10:Singleton=0], human[&&NHX:Genes=9:Identical=9:Duplicated=0:Lost=1:Novel=0:Singleton=0])"human/mouse/mouse2"[&&NHX:Genes=10:Identical=3:Duplicated=0:Lost=0:Novel=7:Singleton=0])LUCA[&&NHX:Genes=3];

I can print an ascii tree like so:

```
/-dog, 3, 0, 0, 0
|
-LUCA /-mouse2, 19, 0, 0, 0
| /mouse_mouse2, 9, 0, 1, 0
\"human/mouse/mouse2", 3, 0, 0, 0 \-mouse, 17, 4, 0, 0
|
\-human, 9, 0, 1, 0
```

Is there a way to do this as a phylogenetic tree image e.g. using the `TreeStyle()` function?

Thanks.

Answer: Yes, each graphical feature associated to a node is called a Node Face. You can add any number of faces to each node, and you can control de position of each face relative to its node. Then call the tree.render(...) function to get an image. There are many type of faces and parameters to customize the image, check this part of the tutorial

A very simple example would like this:

from ete3 import Tree, TextFace, AttrFace
tree = Tree("((A, B), C);")
nodeA = tree & ("A")
nodeAB = tree.get_common_ancestor("A", "B")
nodeA.add_face(TextFace("mouse"), column=1, position="branch-right")
nodeA.add_face(TextFace("/human"), column=2, position="branch-right")
nodeAB.add_face(TextFace("mouse"), column=0, position="branch-top")
nodeAB.add_face(TextFace("/mouse"), column=1, position="branch-top")
nodeAB.add_face(AttrFace("dist", text_prefix="/"), column=2, position="branch-top")
tree.render("tree.png")
Question: Hi all,

We store our alignment files on aws s3. I would like to be able to open them with IGV without needing to download them completely, but I can't find an optimal solution.

If I get a pre-signed url it works but it's not convenient.

I try to follow [this][1] but I can't do it at all. I think the problem comes from the configuration file at the end, I don't know exactly what is needed or not. Knowing that we do not have a connection by google, I would just like to use our aws identifiers. When i click on login on igv and add my longin and password i receive a "succes" answer; but then nothing happen on igv, the menu "amazon --> Load from s3 bucket is still inaccessible.

i use this json (i replace the XXXX by the corresponding value):

{
"apiKey": "",

"project_id": "igv",

"auth_provider": "Amazon",
"aws_region": "us-east-1",
"scope": "email",
"redirect_uris": [
"http://localhost:60151/oauthCallback"
],

"client_id": "XXXXXX",
"client_secret": "XXXXX",
"authorization_endpoint": "XXXXXX",
"token_endpoint": "XXXX",
"aws_cognito_pool_id": "XXXX",
"aws_cognito_fed_pool_id": "",
"aws_cognito_role_arn": "XXXXX"
}

[1]: https://umccr.org/blog/igv-amazon-backend-setup/

Answer: if people come across this post here is some new information:

The last version of igv (2.13.0) add a native support for s3. You just have to get your credentials (programmatic acces) acces key and secret Key, put your credentials files in the good directory (usually it's ~/.aws/credentials ) and now you can click on the Amazon tab in IGV to load your files.
Question: Dear all,

I want to find differential expressed miRNAs between tumor and normal tissues in TCGA data sets.

However when I download whole miRNA transcriptomic data sets, the number of samples belong to adjacent normal tissues or healthy individuals was extremely low(52 tumor samples against only 3 normal samples for Asians, for example).

I want to know if I have downloaded the miRNA data sets correctly.

Looking forward your advice

Nazanin

Answer: Yes, you have most likely downloaded them correctly. For RNA-seq and micro-RNA-seq, the number of normals that were sequenced was minimal. For some of the TCGA cancers, there were no (zero) normals profiled for gene expression. One thing to remember is that most of the 'normal' samples are just from the tissue surrounding the primary tumour, so, they are not representative of true normal tissue.

The TCGA dataset that has most matched Tumour-Normal is breast cancer (BRCA), with ~108 pairings.

You may have to look for other public data if you want to increase your number of normals.

Kevin
Question: Hi everyone. I´m on cytoscape with an undirected network and I want to plot centralities (degree/betweenness), but the values of Betweenness centrality and betweenness unDir are really divergent, and I don´t know which value should I use.

Answer: When calculating betweenness centrality and several other network metrics, it matters if the network is directed or not directed. Betweenness centrality of nodes can be calculated in Cytoscape in at least two different ways: by using the built-in NetworkAnalyzer (in the tools menu) or by using the CentiScaPe app. The reason why this matters is that betweenness comes down to counting how many (or which fraction of) paths go through a given node. For a directed network, the paths must follow the correct directions of the edges. For an undirected network, all edges are considered bidirectional and paths can thus go either way through an edge.

NetworkAnalyzer will prompt the user whether the network should be considered directed or not directed and calculate metrics accordingly, putting the betweenness centrality of each node in a column called *BetweennessCentrality*. In CentiScaPe the user also has to select, but the resulting columns have different names (*Betweenness Dir* and *Betweeness unDir* for directed and undirected networks, respectively).

If your network is undirected, you should always use the undirected version of the metrics. The directed version will be meaningless since it involves assigning arbitrary directions to undirected edges. In case you wonder, the edges will always be considered as pointing from the column you specified as *Source* to the one you specified as *Target* when importing the network.
Question: I was wondering if it is possible to put a title to a `DimPlot()` which is splitted into two by `group.by`.

in my example I am using this command to plot my data

```
DimPlot(seu.Merged, reduction = "tsne", label = TRUE, group.by=c("orig.ident", "seurat_clusters")) + ggtitle('This is a very long title, which should nbe spread over both plots, grouped by Ident (l) and clusters (r)')
```

But the plot looks like that:

![enter image description here][1]

Is there a way to spread the title over both parts of the plot?

[1]: https://i.postimg.cc/fLJ9HQmG/Screen-Shot-2020-09-24-at-7-54-05-AM.png

Answer: *DimPlot* is calling *patchwork* to bind two plots together, so it is better to use *patchwork* to add title, too, try below:

myPlot + plot_annotation(title = 'very very long title, very very long title')

See patchwork manuals for more info:

- [patchwork: Titles, subtitles and captions](https://patchwork.data-imaginist.com/articles/guides/annotation.html#titles-subtitles-and-captions)
Question: I'm trying to run samtools sort in a folder containing about 380 bam files.

I've tried the following command and nothing happens. No errors, it just returns blank. Can anyone help me with what I've done wrong? Thank you :)

samtools sort -@ 7 *.bam > *.sorted.bam

Answer: Try this;

for f in *.bam; do filename="${f%%.*}"; samtools sort -@ 7 $f ${filename}.sorted.bam; done

Or

for f in *.bam; do samtools sort -@ 7 $f ${f/.bam/sorted.bam}; done
Question: I have 25 samples from the TCGA, which contain RNA sequencing expression data from 25 different clinical cancer tumor biopsies. I want to cluster them based on similar expression. The problem is that there are no conditions or replicates to build an experiment design to feed into DESeq or edgeR. I also tried things like perls KMeans library, the R built in kmeans() etc. The problem is that for each sample I have 25K expresssion values (25K genes), so the feature vectors are very large and I don't think I am getting anything that is useful.

Does anyone have any advice on clustering data sets with very large feature vectors and/or clustering expression data without biological conditions.

Thanks,

Kyle

Answer: I'd suggest - first, have a hypothesis. Given what you know about the TCGA biopsies and where they came from, what's your expectation of what a clustering would look like? To get a sense of the data, first make sure the counts are suitably normalized. You could filter the genes to a smaller subset with the highest variance across samples on the log-scale and start clustering with a small set of those genes. Perhaps a few hundred or a thousand, something that is easy to visualize. I would start with using correlation as the similarity measure. Do you you know if the samples came from different clinical sites? Are those sites reflected in the initial clusters? Or the biopsy tissue source? Once you have an initial picture based on a subset of the most variable genes, depending on what you find, you may wish to expand out to include more of the genes, to see what new clusters emerge, if any. Be aware that low-count mRNAs may contribute more noise than signal to your clustering.
Question: Hello there

I am running samtools mpileup and get an output for the requested positions. My problem is that if there are multiple alternative bases, it's not possible to extract the alternative base count for each.

The command line that I use:

samtools mpileup -u -Q 0 -t DP4 -t DV -t DP -f GRCh37.fa --positions positions.bed SAMPLE_1-ready.bam SAMPLE_2-ready.bam SAMPLE_3-ready.bam | bcftools view > SAMPLE-variant-sites.vcf

An example output:

12	25398284	.	C	G,T,	0	.	DP=612;I16=351,184,53,24,19737,1.05346e+06,2839,140417,32100,1.926e+06,4620,277200,9893,221105,1557,36041;QS=2.58612,0.40494,0.0089373,0;VDB=0.0247352;SGB=4.66664;RPB=0.97648;MQB=1;MQSB=1;BQB=0.965306;MQ0F=0	PL:DP:DV:DP4	0,255,255,255,255,255,255,255,255,255:208:3:133,72,2,1	255,0,255,255,255,255,255,255,255,255:186:74:74,38,51,23	0,255,255,255,255,255,255,255,255,255:218:0:144,74,0,0

As you may see, there are G as well as T alternative reads . And if we look at the data from SAMPLE_1, we can see that there is only one DP4 field: `133,72,2,1`. So we know that three bases are either G or T in that position but we do not know **how bases are G** and **how many are T**

How would I achieve this??

Answer: You'll need to update to the latest samtools/bcftools (1.3) and then there are new mpileup tag (-t) options -- AD, ADF, ADR which will give per-allele depth (AD) per-allele depth on the forward (ADF) and reverse (ADR) strands.