cannin/biostars_qa huggingface dataset:
Question: Hi,

I want to download some publicly available differential gene expression datasets for my work. Not sure where to start with and it will be great to get some tips on where to look for.

I checked the GEO for datasets and found files in soft format, but I am looking for files in tab-delimited ASCII text format (similar to the ones used for GSEA)

Answer: Downloading the differential expression data (AKA the actual research finding of a paper) usually turns out to be surprisingly convoluted. It shows the sorry state of the bioinformatics data distribution.

The SOFT files are in fact, ASCII text and tabular format, but you will need to do some parsing to make sense of (drop the initial header lines at the very least). That being said SOFT files represent measures and will not contain other all information that is typically used in a differential expression study.

Data tables with fold changes, p-values, corrected p-values etc are usually only available in the supplementary information published with the paper, usually will be non-standardized, interpreting them requires some manual intervention.
Question: Hi all,

I am not understanding how to download a full matrix of a TCGA subset (e.g. Colon Adenocarcinoma).

I selected Download Data > Data matrix > COAD Data matrix but I am getting a lot of files (level 3).

I want only the matrix of genes x samples. Have I to assemble it or can I download it somewhere?

Thanks

Answer: Hi,

you were searching at the right place. Indeed you will find plenty of txt files when you select "RNA-Seq" in the data matrix. There are gene expression files for each individual seperately. In your case, i assume you want the files endlich with "gene.quantification.txt" as the contain gene and RPKM.

Another possibility would be https://genome-cancer.ucsc.edu/proj/site/hgHeatmap/. Select your tumor and analysis type of interest (Buton "Add datasets") and on the left side (mouse over) you will find a symbol which downloads the data.

I hope that helps,

Sebastian
Question: I want to extract **gene name** , **gene start position** and **gene stop position** from the fasta header of the fasta file. I have tried to extract based on the position but those locations are not consistent. Is there any other way to extract them ?

This is what I have tried so far.

#I have a vector of these file names. Here I have just one element

names1 =>"lcl|NC_005336.1_cds_NP_957781.1_1 [locus_tag=ORFVgORF001] [db_xref=GeneID:2947687]
[protein=ORF001 hypothetical protein] [protein_id=NP_957781.1] [location=complement(3162..3611)] [gbkey=CDS]"

#Then I extracted words from the string list

string_list1 <- str_extract_all(names1, boundary("word"))

#result
string_list1[1]

[[1]]
[1] "lcl" "NC_005336.1_cds_NP_957781.1_1"
[3] "locus_tag" "ORFVgORF001"
[5] "db_xref" "GeneID"
[7] "2947687" "protein"
[9] "ORF001" "hypothetical"
[11] "protein" "protein_id"
[13] "NP_957781.1" "location"
[15] "complement" "3162"
[17] "3611" "gbkey"
[19] "CDS"

So, I was trying to extract 4th ,16th and 17th element from this list. It works for this particular example. This does not work for other headers where these positions are different. Usually, gene name is consistently present at the 4th position. But, the start and stop location differ among the fasta headers. So, this strategy is not working and I can't think of any other strategy.

Answer: Here is the start:

# example data
x <- c("lcl|NC_005336.1_cds_NP_957781.1_1 [locus_tag=ORFVgORF001] [db_xref=GeneID:2947687] [protein=ORF001 hypothetical protein] [protein_id=NP_957781.1] [location=complement(3162..3611)] [gbkey=CDS]",
"lcl|NC_001111_NP_999_1 [locus_tag=Test001] [db_xref=GeneID:2947687] [protein=ORF001 hypothetical protein] [protein_id=NP_957781.1] [gbkey=CDS]")

f1 <- function(x, pattern){
lapply(strsplit(x, " "), function(i){
grep(pattern, i, value = TRUE)
})
}

f1(x, "locus_tag")
# [[1]]
# [1] "[locus_tag=ORFVgORF001]"
#
# [[2]]
# [1] "[locus_tag=Test001]"
f1(x, "location")
# [[1]]
# [1] "[location=complement(3162..3611)]"
#
# [[2]]
# character(0)
Question: Dear all, greetings

i'd like to ask you for a piece of advise please : we have 3 scRNA-seq samples that were sequenced at different depths (200 mil reads, or 800 mil reads, 900 mil reads), and consequently, we do see :

-- distinct numbers of cells, and

-- (on average) distinct number of genes/cell, depending on the sample

would the integration of these samples with CELLRANGER AGGR be a good approach (it does normalize the samples too), followed by standard analysis of the AGGREGATED SAMPLES with SEURAT, or SimpleSingleCell pipeline ?

thank you very much,

-- bogdan

Answer: I would:

**Option 1:**
> Independently quantify genes of each sample. --> Normalize to 10,000
> reads per cell (Default in most scRNA analysis) / **ScTransform** --> Transform the
> matrix to square root (instead of log2(counts+1) --> Merge the three
> matrices (`cbind`)--> Remove genes that are lowly expressed in less
> than 1 or 2 or 5% of the cells --> Use combat to remove batch effects
> (here three batches) --> import the matrix to Seurat --> Skip
> Normalisation --> PCA/UMAPClustering etc... I am pretty sure the cells
> will be clustered by cell-types rather than samples.

If you want to see gene expression changes across clusters, I would introduce an extra step of imputation. So It would be:

> Independently quantify genes of each sample. --> Normalize to 10,000
> reads per cell (Default in most scRNA analysis) --> Transform the
> matrix to square root (instead of log2(counts+1) / **ScTransform**--> Merge the three
> matrices --> Remove genes that are lowly expressed in less than 1 or 2
> or 5% of the cells --> Use combat to remove batch effects (here three
> batches) --> `Impute gene expression` (For example [MAGIC][1]) -->
> import the matrix to Seurat --> Skip Normalisation -->
> PCA/UMAPClustering etc...

Take average of gene expression for each cluster and calculate a cluster specificity score ([Tau Score for example][2]) and them take genes with Tau score more than 0.5 or 0.3 and Perform K-means clustering of averaged gene expression across clusters to pick markers.

**Option 2:**
Use Seurat (v3) CCA analysis to [integrate datasets][3]. Straightforward. It performs **SCTransform** instead of Library size normalisation, which seems to be better for scRNA data but it depends on end goal.

**Option 3:**
Or if you want to use seurat default differential analysis, start with raw counts but use the knn graph from above analysis and proceed with typical marker analysis or differential gene expression analysis.

Its all custom analysis but works pretty well and its fun.

[1]: https://www.cell.com/cell/fulltext/S0092-8674(18)30724-4?_returnURL=https%3A%2F%2Flinkinghub.elsevier.com%2Fretrieve%2Fpii%2FS0092867418307244%3Fshowall%3Dtrue
[2]: https://academic.oup.com/bib/article/18/2/205/2562739
[3]: https://satijalab.org/seurat/v3.1/integration.html
Question: Hi,

I would like to know how could I download all the introns (in FASTA) of a species from Ensembl via web.

Answer: Hi,

As far as I know, we don't store the intron sequences explicitly anywhere.

Here's a piece of Perl code that uses the Ensembl Perl API to fetch all intron sequences for the transcripts overlapping a particular region on the first human chromosome. It can easily be modified to fetch all transcripts in the species and to dump the sequence to a file instead of to the screen:

use strict;
use warnings;

use Bio::EnsEMBL::Registry;
use Bio::EnsEMBL::Utils::SeqDumper;

my $registry = 'Bio::EnsEMBL::Registry';

$registry->load_registry_from_db( '-host' => 'ensembldb.ensembl.org',
'-port' => '5306',
'-user' => 'anonymous',
'-db_version' => '63' );

my $sa = $registry->get_adaptor( 'Human', 'Core', 'Slice' );
my $slice = $sa->fetch_by_region( 'Chromosome', '1', 12_000, 13_000 );

my $dumper = Bio::EnsEMBL::Utils::SeqDumper->new();

foreach my $transcript ( @{ $slice->get_all_Transcripts() } ) {
foreach my $intron ( @{ $transcript->get_all_Introns() } ) {
$dumper->dump( $intron->feature_Slice(), 'FASTA' );
}
}

I hope this helps.
Question: Hi all,

I'm analysing some microarray data and just calculated p-values using a t-test between 2 conditions (each with 3 replicates; total number of genes 14000). My p-values and adjusted p-values are coming out exactly the same. I've tried both FDR and bonferroni for correction using the p.adjust function in R and the results are the same. Is it possible to have the exact same p-values and adjusted p-values. I'm getting a bit doubtful. Any thoughts.

Thanks!!!

Answer: It's impossible to have the same p-value and adjusted p-value after bonferroni correction with 14000 genes. There must be an error with the code at some point, please post it.

Also, you want to use limma rather than directly doing a T-test.
Question: Hi,

I have two lists of genes

mycounts <- read.csv("geneID and length.csv", header = T, sep = "\t", stringsAsFactors = FALSE)

```
> colnames(mycounts)
[1] "genesID.geneslength"

> head(mycounts[1:4,])
[1] "R0010W,1272" "R0020C,1122" "R0030W,546" "R0040C,891"

> dim(mycounts)
[1] 7130 1
```

mycounts1 <- read.table("read.txt", header = T, sep = "\t", stringsAsFactors = FALSE)

```
> dim(mycounts1)
[1] 5961 1
> colnames(mycounts1)
[1] "Freq"
```

How I can have only genes in my read file in my genes file? I mean genes file has 7130 that I only need 5961 of them

May you help me please? Thank you

Answer: See `%in%` operator, e.g.:

data1_subset <- data1[ data1$genes %in% data2$genes, ]

Or we can use `merge`, e.g.:

data1_merge <- merge(data1, data2, by = "gene")
Question: I have a long list of species names.

I was wondering whether a tool was available online, which could fetch the NCBI IDs?

Answer: try [taxonkit](https://github.com/shenwei356/taxonkit), providing single binary files for Windows/Linux. It's really fast.

$ cat names.txt
Homo sapiens
Akkermansia muciniphila ATCC BAA-835
Akkermansia muciniphila
Mouse Intracisternal A-particle
Wei Shen
uncultured murine large bowel bacterium BAC 54B
Croceibacter phage P2559Y

$ time taxonkit name2taxid names.txt
[INFO] parsing names file: /home/shenwei/.taxonkit/names.dmp
[INFO] 1587755 names parsed
Homo sapiens 9606
Akkermansia muciniphila ATCC BAA-835 349741
Akkermansia muciniphila 239935
Mouse Intracisternal A-particle 11932
Wei Shen
uncultured murine large bowel bacterium BAC 54B 314101
Croceibacter phage P2559Y 1327037

real 0m7.217s
user 0m10.682s
sys 0m0.471s

$ time taxonkit name2taxid names.txt | taxonkit lineage -i 2
Homo sapiens 9606 cellular organisms;Eukaryota;Opisthokonta;Metazoa;Eumetazoa;Bilateria;Deuterostomia;Chordata;Craniata;Vertebrata;Gnathostomata;Teleostomi;Euteleostomi;Sarcopterygii;Dipnotetrapodomorpha;Tetrapoda;Amniota;Mammalia;Theria;Eutheria;Boreoeutheria;Euarchontoglires;Primates;Haplorrhini;Simiiformes;Catarrhini;Hominoidea;Hominidae;Homininae;Homo;Homo sapiens
Akkermansia muciniphila ATCC BAA-835 349741 cellular organisms;Bacteria;PVC group;Verrucomicrobia;Verrucomicrobiae;Verrucomicrobiales;Akkermansiaceae;Akkermansia;Akkermansia muciniphila;Akkermansia muciniphila ATCC BAA-835
...

real 0m13.151s
user 0m23.251s
sys 0m0.948s
Question: Does anyone have a simple solution to downloading all the refseq genomes for a particular taxon?

Using `ncbi-genome-download` its possible to specify the species or genus TaxIDs and download them, but apparently you can't go higher up the taxonomic ranks (even though enterobacteria has a TaxID of 543 for instance).

If anyone knows of a way to download all the Enterobacteria I'm all ears.

Alternatively, if there is a method of extracting the species TaxIDs from the Enterobacterial taxid in NCBI such that I can pass them all directly to `ncbi-genome-download` that would work too.

Answer: Maybe this would do the trick!

esearch -db genome -query "txid543 [Organism]"|elink -target nuccore|efilter -query "RefSeq"|efetch -format fasta
Question: Hello all,

I wanted to get a poll of what distros are currently popular with researchers working in bioinformatics and related biological data mining fields.

Mainly:

What linux distribution do you currently (or have previously) use(d)?

What software and language (packs) do you use on a daily basis? -

Where does this distro outperform the rest? Where does it fall short?

I'll start; I currently use Ubuntu and have for the past two years. I mainly use vim with BioPerl scripts. Ubuntu, which is debian based makes installing foreign libraries almost painless. It falls short when it comes to software updates.

Answer: BioLinux works very well. It is equipped with a lot of bioinformatics-related software, and is based on a Ubuntu system. You can also use an Ubuntu system directly, and use BioLinux as a repository for bioinformatics software.
Question: Is there a Picard Tool to filter one or more reads through the following :

FLAG (second column)

Mapping Quality (MAPQ 5th Column)

TLEN (9th Column)

??

Answer: You can use BamTools with its filter utility. You need to write a .json file like this one:

{
"isMapped" : "true",
"mapQuality" : ">50",
"isPaired" : "true"
}

And run command

bamtools filter -in in.bam -out out.bam -script filter.json

See a tutorial here.

I don't see a way to extract TLEN with it, though.
Question: Sorry if this might be a trivial question!

I read a lot about this until I got lost. I need to download wgs VCF file from the 1000 genomes ftp site. I need the snps (snvs and indels), most importantly, I need to have the individual genotypes of all the persons involved.

so for example, this file :

[ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.wgs.phase3_shapeit2_mvncall_integrated_v5b.20130502.sites.vcf.gz][1]

which was referenced many times on biostars, does not contain individual genotypes. I need something similar to what [those files][2] contain. Is there one global file containing snps/indels for wgs data including genotypes of the various samples ?

thanks!

[1]: http://%20ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.wgs.phase3_shapeit2_mvncall_integrated_v5b.20130502.sites.vcf.gz
[2]: ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/phase1/analysis_results/integrated_call_sets/

Answer: You can download the entire data per chromosome (chr1-22 & chrX) —including individual genotypes for both indels and SNPs— using this code:

prefix="ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr" ;

suffix=".phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz" ;

for chr in {1..22} X; do
wget $prefix$chr$suffix $prefix$chr$suffix.tbi ;
done

From: https://www.biostars.org/p/271694/

Kevin
Question: I am an IT guy developing a database for genotype data of samples and I am confused when I look up a snp in NCBI, the genes listed for that SNP sometimes have multiple consequences (intron, missense, upstream, etc.) but on the data that I received only one of those consequences is listed.

Answer: Hello ,

you have to get familiar with the concept of [transcripts][1] to understand this. In short: The sequence of the **gene** can encode for multiple (slightly) different **protein** sequences.

fin swimmer

[1]: https://www.biostars.org/p/244850/
Question: Hi guys,

First I just want to say, I know this has been asked numerous times and in a number of places. However my confusion has increased progressively!

If someone could set the record straight, please, on how to generate a TMM normalized counts using EdgeR, I would be incredibly grateful!

This post has two answers that look like disagree with one another... https://www.biostars.org/p/317701/

The main reason I want to generate TMM matrix from my count matrix is to compare gene expression levels within samples (see the level of multiple cell-type markers within a single sample) in a bar plot, and then generate the same bar plot for every other respective sample to compare between samples.

I think using TMM normalized counts will allow me to compare between samples and within samples according to here: https://hbctraining.github.io/Training-modules/planning_successful_rnaseq/lessons/sample_level_QC.html

I plan to use something like this to generate the bar plots: https://www.biostars.org/p/340485/#341008

Any help would be appreciated.

Thank you in advance :)

Answer: The default normalization in edgeR can be broken down to two steps:

----------

1) normalization by library size. That is simply the correction for
read depth. While this may probably be good enough when there are no
widespread changes in library compisition (=samples are very similar
and only very few genes are differential), this often is not good
enough. See for an example my answer here (https://www.biostars.org/p/9465851/#9465854) using GTEx data where I
compare pancreas and lung transcriptomes, so one would expect
**notably** different gene expression profiles. As you'll see plain per-million scaling results in biased normalized counts while TMM manages to properly center the bulk of genes at y=0 in the MA-plot.

----------

2) the introduction of normalization factors that correct the library
size-scaled values for the compositional component. This here is
what the Trimmed Mean of M-values (TMM) does. For technical details
see the original paper by [Robinson & Oshlack in Genome Biology from
2010][1].

----------

Points 1) and 2) are then combined to calculate the **effective library size** which is then used to divide the raw counts by to obtain normalized counts, also often referred to as TMM-normalized counts or `cpm`.

In practice:

#/ make the DGEList:
y <- DGEList(...)

#/ calculate TMM normalization factors:
y <- calcNormFactors(y)

#/ get the normalized counts:
cpms <- cpm(y, log=FALSE)

The `cpm` function uses the normalization factors (given that `calcNormFactors` was run on that DGEList) internally.
If not, then `cpm` just return the plain per-million scaled factors.

[1]: https://genomebiology.biomedcentral.com/articles/10.1186/gb-2010-11-3-r25
Question: I am trying to use PAML for the first time, and it's thrown up a few errors I am struggling to fix. The manual is particularly tricky to read, but I've solved most of the other problems I've come across. I am aiming to use it to calculate dN/dS for sets of genes. I used RAxML to generate my trees in Newick format, but the error message I get after executing codeml is this:

> Error: need branch labels in the tree for the model..

My tree looks like this:

(chrI_BS45_CA_R_ENSGACG00000022906:2,chrI_BS49_CA_R_ENSGACG00000022906:1,chrI_BS43_CA_R_ENSGACG00000022906:1,chrI_BS47_CA_R_ENSGACG00000022906:1,chrI_BS51_CA_R_ENSGACG00000022906:1,chrI_BS53_CA_R_ENSGACG00000022906:1):1.0;

And my control file looks like this:

seqfile = /data/12_PAML/Input/Chr1/CA_R_ENSGACG00000022906_pruned.fa
outfile = /data/12_PAML/Output/Chr1/CA_R_ENSGACG00000022906_Codeml_Output.txt
treefile = /data/12_PAML/Input/Chr1/CA_R_ENSGACG00000022906_pruned.tre

noisy = 9
verbose = 1
runmode = 0

seqtype = 1
CodonFreq = 0
model = 2
NSsites = 2
icode = 0
fix_kappa = 0
kappa = 1
fix_omega = 0
omega = 1
cleandata = 1

I've tried looking for similar answers and Ziheng has posted about similar problems saying to look at the examples folder. In it you get a trees like this:

((1,2) #1, ((3,4), 5), (6,7) );
or
((_10_H._fulgens, (((__4_H._kamtschatkana, (__1_H._rufescens,
(__2_H._sorenseni, __3_H._walallensis))) #1, (__5_H._sieboldii,
(__6_H._discus_hannai, __7_H._gigantea))), (__8_H._corrugata,
__9_H._cracherodii))), _25_H._iris, ((_17_H._pustulata,
(_24_H._t.coccinea, _23_H._t.tuberculata)), (_22_H._australis,
((_18_H._midae, ((_11_H._roei, (_12_H._scalaris, _13_H._laevigata)),
(_14_H._cyclobates, (_15_H._rubra, _16_H._conicopora)))),
((_19_H._ovina, _21_H._varia), _20_H._diversicolor)))));

I'm struggling to see the difference. I've tried removing the floating point values and running it again, but to no avail.
Any help would be very much appreciated.
And thank you for your time.

Answer: Turns out this was a very simple answer. For using any of the PAML models with selection, you need to designate which branch(es) on your gene tree or phylogeny where selection differs. The notation is #1 or #*n* if there are *n* branches with different strengths of selection.
Question: Hi,

I have mapped (paired-end) RNA-seq reads (so BAM file information) and I want to extract the nucleotide in the genome immediately prior to the (forward) mapping (essentially to test whether a RNAse digest has worked correctly).

As it is in BAM format, then extracting the position of the read (chromosome and location) is straightforward (just samtools view it, and process the line output for the relevant field), but I was wondering if there was a quicker way than brute-force (i.e. extracting, sorting, then reading from the genome FASTA file).

Answer: Using my tool [**samslop**](https://github.com/lindenb/jvarkit/wiki/SamSlop) and [**bioalcidae**](https://github.com/lindenb/jvarkit/wiki/BioAlcidae):

$ java -jar dist/samslop.jar -c -m 1 -M 1 -r ref.fa S1.bam |\
java -jar dist/bioalcidae.jar -F SAM -e 'while(iter.hasNext()) {var read=iter.next(); if(read.getReadUnmappedFlag()) continue; out.println(read.getReferenceName()+":"+read.getAlignmentStart()+"-"+read.getAlignmentEnd()+" "+read.getReadString().substr(0,1)+"/"+read.getReadString().substr(read.getReadLength()-1));}' | uniq | head

rotavirus:1-71 G/A
rotavirus:1-72 G/A
rotavirus:1-53 G/A
rotavirus:1-72 G/A
rotavirus:1-65 G/A
rotavirus:1-68 G/A
rotavirus:2-73 G/T
rotavirus:3-74 C/A
rotavirus:3-69 C/T
rotavirus:3-74 C/A
Question: I calculate the pangenome for a set of bacterial genomes I am working with.
I have a classical pangenome matrix that looks like this:

speciesA speciesB ... speciesZ
gene1 1 2 0
gene2 0 1 0
gene3 1 1 1
gene4 1 1 2

Now I would like to map the gen gain and loss in a phylogenetic tree.
Something like this:

![enter image description here][1]

I guess I could just record all families of genes shared among the different species in each branch of the tree, but I believe there is a more elegant way of doing that.

Any suggestion?
[1]: http://www.nature.com/nature/journal/v510/n7503/images/nature13400-f3.jpg

Answer: I guess you are looking for a software like Count, [that performs ancestral reconstruction of gene family sizes over a given tree][1]. But as [abascalfederico][2] mentioned, you should be very careful with the possibility of HGT.

[1]: http://bioinformatics.oxfordjournals.org/content/26/15/1910.full "doi: 10.1093/bioinformatics/btq315"
[2]: https://www.biostars.org/p/180593/#180615
Question: So here is the result from cuffdiff, I only printed out certain columns for read:

```
test_id gene value_1 value_2 significant
XLOC_226942 - 3.73672 0 yes
XLOC_227336 - 6.80194 0 yes
XLOC_227426 - 0 111172 no
XLOC_227811 - 0 154.721 no
XLOC_228955 - 0 9.80122 yes
XLOC_229214 - 6.61332 0 yes
XLOC_229754 - 8000.99 0 no
XLOC_231416 - 0 714.61 no
XLOC_231570 - 11.8838 0 yes
XLOC_231696 - 109.917 0 no
XLOC_231746 - 0 946.635 no
XLOC_231910 - 2.8676 0 yes
XLOC_231945 - 7.25977 0 yes
XLOC_232519 - 2.42438 0 yes
XLOC_232650 - 126.899 0 no
XLOC_232742 - 25.7937 0 yes
XLOC_233457 - 163.295 0 no
```

So why it shows no significant difference if one of value equal to 0 but another is over 100? Isn't that a huge different? And if one is 0 and another is less than 100, it considered significant. Do you guys have the explanation for this? Do I also need to consider those as significant expressed genes?

Thanks so much!

Answer: This is a known issue with Cuffdiff, and one reason people generally prefer other programs (DESeq and edgeR) instead.