cannin/biostars_qa huggingface dataset:
Question: Hi, I've a doubt for which I've no answer.

The fourth field of the TCGA barcode contains the information about the sample, either primary tumor, metastatic, etc

For instance `TCGA-93-A4JP-01A-11H-A24S-13`, this is a primary tumor, but what does the A next to the 02 means? Sometimes I've sen other letter than A such C, E...

Thanks

Answer: **Edit April 6th, 2020:**

Since I posted this answer, this page seems to have come online: TCGA Barcode

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

# TCGA Barcode

# Code tables

Go here for the code tables for further in depth details:

- https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables

Kevin
Question: Dear all,
I am a rookie in data analysis and stuck with my results dnt know how to interpret them.

I started with 7 metagenomic assemblies of different species of Azolla fern. The aim was to identify bacteria in leaf ecosytem of azolla different species. Out hypotheisis was, if there are similar bacteria which repeat within the azollas different species, they will cluster together when their genomes will be plotted in dendrogram or a tree.

The method used `spades` to get assemblies, `BWA` was used to do backmapping, `samtools` for sorting, `metabat` for binning and checkm for to see completeness and contamination of bins.

Then `prokka` was used to annotate the genomes and uniport ids were obtained and table was made of all uniport id of all the bins. the table was changed to binary table and then used to create a dendrogram in R.

The dendrogram and then tree made by using dendrogram in fig tree. In the tree i observed that the bacteria are clustering according to the metagenomic sample or plant host not on the basis of their similar taxonomical name eg rhizobiales is clustering with burkholderiales of same metagenomic assembly but not with rhizobiales of other host plant assembly.

Im on the dead end how to intrepret these results and what can i deduce from it.
and are there other ways to improve my approach? Can i compare similar taxonomical bins directly of different metagenomic assemblies any suggestions will be valuable.

kind regards

manpy student utrecht university holland

Answer: An alternative approach might be to calculate `mash` distances between all your genomes directly. Assuming your taxonomic assignments are largely correct, this would save you having to mess about with extracting the uniprot IDs, binary matrix etc.

Since you can’t align that many full genomes, mash distances are a good surrogate for genome similarity.

You can then draw your trees using pairwise mash distances among all your genomes.

I’m not 100% sure how much you’ll need to polish your genomes etc. You may need to reorder contigs and perhaps concatenate but if you read up on it I’m sure it’ll become clear.
Question: Hi everyone!

I am trying to find the best way to make 2 boxplot for a specific gene from data found in a row for a subset of columns within data frame x.

x dimensions are 634 by 128 columns

Each row is specific to a gene,

Column 1 has gene name, and I want to say look at gene in row#1

columns 2:48 data I want to include in one boxplot

columns 49:128 data I want to include in another boxplot

data frame looks something like this

```
gene accepted_hits_x1.bam accepted_hits_x1.bam etc....
1 AARS1 -6 0 etc....
```

I also want to be able to see each data point that makes up the boxplot plotted in the plot

I am having a problem:

I am running into the problem where my data (residual from mean ... meaning x value - mean) is a series of positive and negative values and it appears that with this plot it is excluding these negative values...

```r
data <- unlist(subset(datavr, gene =="IGF1R", select=2:128))

news <- data.frame(data=data, factor=c(rep(1,47), rep(2,80)))
news$data <- (log10(as.numeric(news$data)) + 1)

g <- ggplot(data=news, aes(x=as.factor(factor), y=data))

g + geom_boxplot() + geom_point(color="purple", size=3) + xlab("A38-41 A38-5 ") + ylab("log10(Residual from Mean)+1") + ggtitle("IGF1R inside region") + theme(plot.title = element_text(face="bold"))
```

The problem is that it keeps giving me error saying that:

```r
Removed 110 rows containing missing values (geom_point)
```

This could be that these values are negative so taking the log10(value)+1?

Answer: ```r
gene_id <- 1 # consider the first gene
data_1 <- your_dataframe[gene_id,2:48]
data_2 <- your_dataframe[gene_id,49:128]
boxplot(data_1,data_2)
```
Question: Hello,

I am looking for an article or video which explains how the normalization in the limma voom is performed. I looked for some articles but I could not find any which explains it clearly. Can someone help me in this regard?

Answer: A good place to start would be to look back to where limma was originally used, i.e., for microarray analyses and microarray normalisation, and i highly recommend Professor John Quackenbush's great paper on this: http://www.cs.cmu.edu/~zivbj/class05/reading/norm.pdf

Then I would read-up on linear modelling and linear regression, as limma is fundamentally based on linear modelling.

Finally, with that under your belt, I would read the actual work by the 'new' limma authors:

- https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4402510/
- https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4937821/

Part of being a great researcher is being able to read and understand your colleagues' published works, but I can understand that entering a new field can take time when there are literally 1000s of articles already out there on the topic.
Question: Dear all:

I have downloaded Geneious software trial version [https://www.geneious.com/][1], but now it is expired. As a fresh level grad student, I can't afford to pay full version of Geneious software. I aware that Geneious has abundant features to analyze NGS data. My expectation is to find tools that relatively similar to Geneious. My major intent is to carry out **phylogenetic analysis**.Can anyone recommend the better alternative for Geneious tools? Is there any free, open source software tools that can be alternative for Geneious? My colleague is not well comfortable with R environment, so we are seeking tools that require less programming input. Can anyone give possible aid or point out which library we can use?

[1]: https://www.geneious.com/

Answer: You can also have a look at [UGENE][1] as an alternative to Geneious.

[1]: http://ugene.net/
Question: Hi all,

I'm a bioinformatician and I often have to deal with RNAseq differential gene expression analysis projects.

I think I understand well the whole process to get from the raw data to the the normalized read counts but, unfortunately, due to my little statistical background, I'm having trouble dealing with the last step of differential expression.

When it comes to simple pairwise comparison between two conditions I understand the process, but when there is more complex comparisons (timelapses, multiple comparisons, including confounding effect ... ) I'm struggling for chosing the relevant design matrices.

I'm curious if anyone would know good tutorials, online courses, books, or any ressources that would allow me to learn how to get better at that.

Thank for your help,

Answer: I've not actually done it myself, but people say good things about Rafael Irizarry's courses. Particularly see weeks 3 and 4 of the Introduction to Linear Models and Matrix Algebra course [here][1] and [chapter 5][2] of his book.

[1]: http://rafalab.github.io/pages/harvardx.html
[2]: http://genomicsclass.github.io/book/
Question: I have two files with genes

File one (with 40000 genes)

Gene 1
Gene 2
Gene 3
Gene b
Gene f
Gene c
Gene r
Gene z

File two (with 39000 genes)

Gene 1
Gene 3
Gene 2
Gene b

I would like to know if there is a command line (with awk or bash) to extract that lines that exist in the one file and not in the two file

Answer: > I would like to know if there is a command line (with awk or bash) to extract that lines that exist in the one file and not in the two file

use comm : http://man7.org/linux/man-pages/man1/comm.1.html

comm -3 <(sort file1.txt) <(sort file2.txt )
Question: I am analysing protein data using Python programming language and Jupyter notebook.

In the Terminal I have put an alias in a hidden file on the home directory entitled .bash_profile, in order to be able to open pymol directly from the Terminal .

alias pymol=/folder/where/pymol/is/located/Applications/PyMOL.app

In Jupyter notebook, the first two commands are executed:

from ipymol import viewer as pymol
pymol.start()

However, the following command gives an attribute error.

pymol.fetch('4MBS') # Fetch PDB

Is there another way to load a file directly from my hardrive for viewing?

The only workaround I have found is using.

import nglview as nv

view = nv.show_file(/path/to/folder/where/protein/located.protein)
view

I appreciate any help in advance as i would like to use pymol in this way, if possible.

Thanks

Answer: The answer to this question is that a bash shell .pml script needs to be made instead of a python .py script.

The shell script can be shown in a Jupyter Notebook by using RawNBConvert setting, which is found underneath the help menu.

The images generated from the .pml script can be viewed using Pymol and snapshot jpeg or png files can be taken of the molecules.

The ipymol does not work and the nglview is not required.
Question: Hi All,

I know there are a few posts here that raise specific questions about RNA-seq library prep. protocols, but I was curious if there's a comprehensive catalog about exactly what protocols exist and exactly what type of data they produce (i.e. to what do the reads in the final FASTA/FASTQ files correspond). Basically, I'm somewhat confused by all of the different aspects of a protocol and how they compose. For example, if a protocol is stranded or not, the relative orientation of the reads, which strands each read (i.e. \1 and \2) comes from. Are the reads that end up in a FASTA/Q file always reported 5' to 3' with respect to the strand from which they are derived? Are mates pointing toward each other reported with respect to the same strand, opposite strands, or both depending on the protocol? What do people mean when they say reads are 'reversed' --- does this mean reverse complemented with respect to the mate, or that the read is actually reported in the 3' to 5' direction (i.e. reversed but not complemented)?

Basically, I'm curious how all of these different variables interact with each other to produce the read sequences that will be used for downstream analysis. If I want to be able to communicate (to another person, or, perhaps as importantly, to a piece of software) all of the details / constrains about which reads should map to which strands in which orientations --- what is the most parsimonious way to do so? What is the minimum amount of information I need to convey? Is there a standard language / specification for representing this information? I'm sorry to ask such a broad question, but I'm a bit overwhelmed and trying to gain a comprehensive understanding of what, exactly, the reads in a file represent in light of the protocol under which they were prepared.

Thanks!

Answer: The relative orientation of a pair of reads will be the same unless you're using mate-pairs or the like. This is true regardless of whether you have a stranded/directional library or not.

Edit: I should add that when this isn't the case with a standard paired-end library then either (A) an error occurred when the library was made (PCR or otherwise), or (B) the sample doesn't match the reference there (i.e., you have a variant), or (C) it's a case of incorrect mapping, or (D) something really strange is going on (perhaps the nature of the experiment would make this likely).
Question: Here is my attempt:

def remove_too_long_reads(bam):

for read in bam:
if read.end - read.start < 60:
yield read

bam = BedTool(input.bam)
bam = BedTool(remove_too_long_reads(bam))

But when I try to use the bam file afterwards like so:

bam.intersect(...

I get an error:

> : Broken pipe The command was:
>
> bedtools intersect -b data/gencode_annotation/tss.bed -a stdin
>
> Things to check: Error in job intersect_data_tss while creating output
> file data/tss/bam/Exp2_9h_PolII.bam. RuleException: KeyError in line
> 31 of
> /local/home/endrebak/code/programmable_epigenetics/rules/tss/intersect_bam_gencode.rules:
> 32 File
> "/local/home/endrebak/code/programmable_epigenetics/rules/tss/intersect_bam_gencode.rules",
> line 31, in __rule_intersect_data_tss File
> "/local/home/endrebak/anaconda3/lib/python3.5/site-packages/pybedtools/bedtool.py",
> line 773, in decorated File
> "/local/home/endrebak/anaconda3/lib/python3.5/site-packages/pybedtools/bedtool.py",
> line 336, in wrapped File
> "/local/home/endrebak/anaconda3/lib/python3.5/site-packages/pybedtools/helpers.py",
> line 373, in call_bedtools

Answer: In order to maintain compatibility with VCF/BED/GFF/GTF functionality, pybedtools converts BAM/SAM reads into `pybedtools.Interval` objects. When just filtering on length as you're doing here, that's unnecessary overhead. Here's a more efficient way to do it with pysam, which also has the benefit of providing access to all the read details if you need to do more sophisticated filtering. You can then use the filtered bam for downstream operations (e.g., histograms of read counts across intervals):

https://gist.github.com/daler/f44c19510860aeb29e8f75432eec8883
Question: I would like to compare the results from different algorithms for clustering data from an RNAseq experiment. One methodology we used is **WGCNA**, which represents each cluster/module by a module **eigengene**.

What I am now looking for is a way to calculate such eigengenes on arbitrary lists of genes and their expression profiles. I want to use that to reduce the clusters generated by other methods to a set of "representative genes" and compare that to the WGCNA output.

Unfortunately, I'm not clear about how to get from a list of genes and their expression profiles to an eigengene. The functions inside the WGCNA R package are deeply tied into the WGCNA analysis and I can't see how to use them on arbitrary data frames of gene expression data. Any hints would be highly welcome!

Answer: Eigengenes were first defined in this [paper][1]. Singular value decomposition (SVD) is what you're looking for. The eigengenes are the right singular vectors of the SVD of the expression matrix. If X is your data with genes as rows and samples as columns, the SVD of X is X=USV' and the eigengenes are defined as the vectors in V. In R,

eigengenes <- svd(X)$v

EDIT: Fixed link to paper.

[1]: http://www.pnas.org/content/97/18/10101
Question: blastp -query homo.faa -db cow.faa -out homo_blast.csv -outfmt 6 -evalue 0.00001 -max_target_seqs 1 -num_threads 32

I got the below error while performing protein blast with the code mentioned above

> BLAST Database error: No alias or index file found for protein
> database [cow.faa] in search path

When I revised the code by specifying the path, the problem was not resolved. When I typed subject instead of db, the problem was solved, but I'm not sure if it's correct.

Answer: When you use a `-db` switch, it is expected that your database will be indexed. This is done using a `makeblastdb` command, and will create several files that will end in `.p??` where ?? stands for two other letters. When searching against a large database, to gain speed you would want to index it and use a `-db` switch. You most likely have not done this step.

When using a `-subject` switch, the database need not be indexed. The result will be the same as when using the `-db` switch, but the search will be slower and this is not recommended for larger databases. It is typically used for comparing two sequences, or when the target database is relatively small.
Question: I try to understand how to chose the optimized pi_hat parameter for a dataset. In many articles, they chose 0.2 as pi_hat, and everything above that is considered to be cryptic relatedness or duplicates.

I've tested IBD on HapMap, the files I use can be found here: ftp://ftp.ncbi.nlm.nih.gov/hapmap/genotypes/2009-01_phaseIII/plink_format/. I first remove all annotated offspring from HapMap. Then I peform IBD to see if it still finds samples with cryptic relatedness to each other. The steps I peform are the following (in PLINK):

1) LD-prune:

plink --file hapmap --indep-pairwise 50 5 0.2
plink --file hapmap --extract plink.prune.in --recode --out hapmap_pruned

(2) IBD:

plink --file hapmap_pruned --genome --min 0.2

The results shows that many cryptic related samples can be found with a pi_hat of 0.2 as threshold, even if all offspring were initially removed. My question is, is this a normal behavior? Or should one increase the pi_hat? How to find out a "good" pi_hat for a custom dataset?

Answer: This is normal behavior for the HapMap data set. See Stevens et al. and a few of the earlier papers he cites that try to identify un-annotated relationship in HapMap. The Stevens paper is using Cotterman coefficients (K1,K2) which are fraction of the genome shared IBD1, and IBD2. As zx8754 mentioned above, 0.125 is third degree, 0.25 is second degree, and 0.5 is first degree, although in practice these thresholds can be too low or high due to consanguinity, and admixture.
Question: I used this piece of biopython codes written by Eric Normandeau to parse several sequences from a Fasta file. The codes worked nicely. I would like to have my parsed output from this format

>DR179241 similar to UniRef100_A5HIY3 Cluster: Thaumatin-like protein; n=1; ...... ATAATCTTTAGATCAGTCATCAATCTCAACAGTATCGCTTTCAATTCTCTTTCATATTGC ATGGAAGTGTGTAAATACAATTAGGGCATTCATTGAGTTGACTTCATTTAAGCGCT......

exactly like this format:

>DR179241
ATAATCTTTAGATCAGTCATCAATCTCAACAGTATCGCTTTCAATTCTCTTTCATATTGC
ATGGAAGTGTGTAAATACAATTAGGGCATTCATTGAGTTGACTTCATTTAAGCGCT......

Could anyone kindly please light me for example how to modify this piece of biopython codes written by Eric Normandeau in order to suit my purpose?

Thank you very much and have a nice day.

Answer: You could do this in R, using seqinR:

library(seqinr)
x=read.fasta("yourfile.fasta")
write.fasta(x,file="newfile.fasta",names=names(x))
Question: I am stuck with a very simple problem. I want to build a Pearson correlation matrix for my microarray dataset. My .cvs file consists of normalized, log-transformed expression values of 18k genes across 36 samples. I want to find the gene-gene Pearson correlation from this matrix using R package. After that, I want to transform the matrix to the form of an edge-list with genes in the first two columns and the value of the correlation in the last column. I was trying out the cor() function in R, but I guess there is some issue with numeric/character values because of which it gives me the error `x has to be numeric`. Kindly give some suggestions as to what way I can read in the file and transform the matrix.

Thanks

Gene sample1 sample2 sample3
A 10 50 78
B 50 45 55
C 70 56 44

Answer: Use my script taxo_bivariate_plot.R. It uses `cor()` as suggested by Phil S. Usage information is as follows:

$ Rscript taxo_bivariate_plot.R --help
Usage: taxo_bivariate_plot.R [options] file
Options:
--ifile=IFILE
CSV file
--opath=OPATH
Output path
--fsize=FSIZE
Font size [default 1.2]
--width=WIDTH
Width of jpeg files [default 800]
--height=HEIGHT
Height of jpeg files [default 800]
--correlation=CORRELATION
Correlation to use: 1=pearson, 2=spearman, 3=kendall [default 1]
--rmode
Mode: TRUE=R mode, FALSE=Q mode [default FALSE]
-h, --help
Show this help message and exit

This script generates bivariate plots with histograms on the diagonals, scatter plots with smooth curves below the diagonals and correlations with significance levels above diagonals. Data file has the following organization:

Var_1 Var_2 Var_3 .. Var_R
Sample_1
Sample_2
Sample_3
...
Sample_N

For example,

$head ENV_pitlatrine.csv
Samples,pH,Temp,TS,VS,VFA,CODt,CODs,perCODsbyt,NH4,Prot,Carbo
T_2_1,7.82,25.1,14.53,71.33,71,874,311,36,3.3,35.4,22
T_2_10,9.08,24.2,37.76,31.52,2,102,9,9,1.2,18.4,43
T_2_12,8.84,25.1,71.11,5.94,1,35,4,10,0.5,0,17
T_2_2,6.49,29.6,13.91,64.93,3.7,389,180,46,6.2,29.3,25
T_2_3,6.46,27.9,29.45,26.85,27.5,161,35,22,2.4,19.4,31
T_2_6,7.69,28.7,65.52,7.03,1.5,57,3,6,0.8,0,14
T_2_7,7.48,29.8,36.03,34.11,1.1,107,9,8,0.7,14.1,28
T_2_9,7.6,25,46.87,19.57,1.1,62,8,13,0.9,7.6,28
T_3_2,7.55,28.8,12.65,51.75,30.9,384,57,15,21.6,33.1,47
$ Rscript taxo_bivariate_plot.R --ifile=ENV_pitlatrine.csv

Will generate the following image:

This Rscript is back-end of my TAXAenv website. You can use `--rmode` to transpose your matrix and change the script accordingly to meet your needs.

Best Wishes,
Umer
Question: I know I might be missing something obvious, but how do I get the accuracy/quality of my imputation after using IMPUTE2 (or any imputation tool)? How is the quality of an imputation measured?

Thank you!

Answer: I actually figured it out. Incase anyone ever has the same question. IMPUTE2 outputs a file called info where it has imputation quality information. Depending on the command line arguments you choose you can get more specific statistics.
https://mathgen.stats.ox.ac.uk/impute/output_file_options.html#info_metric_details
Question: Hello everyone,

I am new in 1000 genomes project data. I want to download all bam files belonging to phase3, can anyone guide me how can I download all of them (from the command line?). Do you have any estimation how long it is going to take?

I want to compute the depth of coverage only for some specific intervals, not the entire genome. Is there any way that I could do it without downloading the data?
I could find this, but not sure if it is relevant to what I want to do?

samtools view -b ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data/HG01375/alignment/HG01375.mapped.ILLUMINA.bwa.CLM.low_coverage.20120522.bam 2:1,000,000-2,000,000 | genomeCoverageBed -ibam stdin -bg > coverage.bg

I would appreciate if anyone could guide me.

Answer: you wrote:

samtools view -b ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data/HG01375/alignment/HG01375.mapped.ILLUMINA.bwa.CLM.low_coverage.20120522.bam 2:1,000,000-2,000,000 | (...)

you want:

samtools -bu view 'http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01375/alignment/HG01375.mapped.ILLUMINA.bwa.CLM.low_coverage.20120522.bam' "2:1000000-2000000" | (...)
Question: Hi,

I have this list of proteins from a new genome project so its pretty much unannotated. However, it's closely related to *C. elegans* so I was thinking of trying to identify the closest *C. elegans* homologues.

What I've been doing right now is doing a protein blast in ncbi with the protein sequences and then taking the top *C. elegans* hit, however, there are far too many sequences to be able to do this one at a time, so I was wondering if there's a way to do it faster/automated/program that does it for me.

Thanks!

Answer: Try using [HMMER][1].

The manual is available [here][2].

IN BRIEF: for each protein sequence in *C. elegans* you make a HMM using `hmmbuild` command. Concatenate all HMM models into a single file to make a database file. You have to use `hmpress` to create additional files in order to search your database. Now you can use either `phmmer` if you want to scan against the database you have just created or `hmmsearch` to scan individual models against the sequences you have. The documentation describes very well what commands you need, but note the subtle differences of scanning model vs set of sequences and set of sequences vs db of models.

If you have access to a parallel environment such as MPI (OpenMPI can usually be installed even on the local machines to take full advantage of multiple cores) then you can build the HMMER with MPI support to increase throughput.

A rough idea of a time in our use was: Building and pressing a database of ~10k models takes 10 mins (ish) scanning a coding sequence against a database of ~10k models takes 2-3 seconds. This is very rough guide that we have used it, which undoubtedly will differ from your use case.

[1]: http://hmmer.janelia.org/
[2]: ftp://selab.janelia.org/pub/software/hmmer3/3.1b1/Userguide.pdf