cannin/biostars_qa huggingface dataset:
Question: Hello Biostar,

I'm performing some Gene Ontology enrichment analysis using R/Bioconductor and, as far as I understand, there are two packages that allow to do this: topGO and GOstats (I'm not interested in packages that compare multiple datasets in this scenario).

The main difference seems to be the statistics used: GOstats uses the hypergeometric distribution (basically the standard way to test for overrepresentation), while topGO allows the user to use a much wider and complete range of algorithms and statistical tests to check for enrichment.

With regards to graphic capabilities, topGO seems to be producing much better graphs than GOstats, but I have only looked at the packages' vignettes.

Has anybody here experience with both? Which one did you end up using? Why?

Also, on a related note, what would be the topGO algorithm/statistics test combination that emulates GOstats' hypergeometrical test? Would it be the classic/fisher combination (used as the most basic example in topGO's vignette)?

Thank you.

Answer: To avoid the following situation:

I tried both GOstats and topGO briefly before switching to topGO for all of my GO enrichment analysis in R. Its been a while but from what I remember it seemed more versatile, i.e. one can employ various default algorithm/test combinations (along with the capability of customized tests).

It seems reasonable that the classic/fisher would emulate GOstats hypergeometric test, though I have not tested this personally.

I found this was a good resource for working with topGO.
Question: Hi,
I am analyzing genomic (RNA-seq) data from Patient Derived Xenograft tumor samples where cancer patient tumors are transplanted and grown in a mouse, harvested, and then extracted DNA and RNA is sequenced. I have never done this before and wondering if human or mouse reference genome should be used. I would guess the human reference genome should be used for alignment but is it possible that there would be some mixed in mouse cells?
Thanks,
- Pankaj

Answer: First, use some wet lab technology to filter out mouse cells prior to DNA/RNA extraction (just google "mouse cell depletion"). Then, in Bioinformatics analysis, map all reads to human and mouse genome/transcriptome in parallel and independently. Reads that map better to mouse are discarded from the analysis and only "human reads" are kept.

As poisonAlien mentioned, some tools exist to remove contamination but you have to carefully evaluate them. For example, some do this based on k-mer counting and can be too sensitive or too specific. We found, that the above described approach gives best results even though it is computationally demanding. In particular, if you want to call variants, this approach removes fewest reads but gives high sensitivity.

My company is specialized in xenograft research, so feel free to contact me for more information.
Question: I have 2 sequences with gaps:

15888 CATG----ACAGAGCGACCCGCG--CACGTTACAAACACTACGCGGGGTGGCCCCGG
42690 CATGCCCGACAGAGCGACCCGCGAACACGTTACAAACACTACG---GGTGGCCCCGG

I would like to encode indels in binary form. Gap is marked with 0 and no gap is 1. I want to receive something like this:

15888 0 0 1
42690 1 1 0

Thanks for any answer.

Answer: As suggested, consider the *the larger question*. But, based on assumptions from your example, looks like an interesting problem.

Here is using *R*, not going to explain the code in detail. In short, find position of `"-"` in all strings, create groups based on consecutive positions, i.e.: make *"indels"* ranges, then compare those position for each string.

foo <- function(x){
xSplit <- strsplit(x, "")

# indels
indels <- unique(
unlist(
lapply(xSplit, function(i){
ix <- which(i == "-")
g <- c(0, cumsum(diff(ix) > 1))
split(ix, g)
}), recursive = FALSE))
# sort by position
indels <- indels[order(sapply(indels, `[`, 1), lengths(indels))]
names(indels) <- paste0("indel_", seq(length(indels)))

# return binary
t(
sapply(xSplit, function(i){
sapply(indels, function(j){
!all(i[ j ] == "-")
})
}) * 1
)
}
# example data
x1 <- c("CATG----ACAGAGCGACCCGCG--CACGTTACAAACACTACGCGGGGTGGCCCCGG",
"CATGCCCGACAGAGCGACCCGCGAACACGTTACAAACACTACG---GGTGGCCCCGG")
x2 <- c("CATG----ACAGAGCGACCCGCG--CACGTTACAAACACTACGCGGGGTGGCCCCGG",
"CATGCCCGACAGAGCGACCCGCGAACACGTTACAAACACTACG---GGTGGCCCCGG",
"CATG----ACAGA---ACCCGCG--CACGTTACAAACACTACGCGGGGTGGCCCCGG")

# test
foo(x1)
# indel_1 indel_2 indel_3
# [1,] 0 0 1
# [2,] 1 1 0

foo(x2)
# indel_1 indel_2 indel_3 indel_4
# [1,] 0 1 0 1
# [2,] 1 1 1 0
# [3,] 0 0 0 1
Question: Hi,

I need to analyze down-sampled data of couple of Rna_Seq full data set (samples,paired-end,fastq).sub-sampling method should work the same for all samples. In the end I compare for example how 5% of full data differs from 10% of full and 20 and 40% as well. (a sample is :ERR188044)
The final graph will depict how amount of data affects the result.

The question is :
How to download the data in these four forms ? shall I firstly download the full data and then downsample ? Or I can directly download down-sampled data.
how to sub-sample data to get only a few number of specific chromosomes ?
how to sub-sample data to get only a percent of whole paired-end reads?

What do you suggest me to do ?

Your advice is appreciated.

Thanks.

Answer: FYI, the term you're looking for to describe this is a "rarefaction curve".

Yes, you will need to download the full dataset and then subsample it (typically a few times) for each percentage you want to plot on the X axis. You can use things like `seqtk` or even `samtools view` to subsample files. Depending on what you need to do, aligning the whole dataset once and then subsampling that will likely turn out to be the quickest strategy.
Question: cd-hit has given me an output that looks just like the example below.

```
>Cluster 0
0	496aa, >SRR5892231.2396932... *
>Cluster 1
0	496aa, >SRR5892231.3763255... *
1	390aa, >SRR5892231.1558909... at 91.03%
>Cluster 2
0	496aa, >SRR5892231.1710795... *
>Cluster 3
0	496aa, >SRR5892231.2083014... *
1	464aa, >SRR5892231.14158... at 91.59%
2	423aa, >SRR5892231.1116524... at 94.56%
3	314aa, >SRR5892231.1717279... at 95.86%
4	268aa, >SRR5892231.2309241... at 99.63%
5	371aa, >SRR5892231.480233... at 99.46%
>Cluster 4
0	496aa, >SRR5892231.3954388... *
1	319aa, >SRR5892231.1752373... at 99.69%
>Cluster 5
0	496aa, >SRR5892231.14746... *
>Cluster 6
0	496aa, >SRR5892231.2340653... *
1	407aa, >SRR5892231.2608197... at 100.00%
2	340aa, >SRR5892231.1216749... at 100.00%
3	345aa, >SRR5892231.3205930... at 92.46%
```
Each line that starts with a > shows the cluster label.
The next line that ends with a ... * shows the representative sequence for the particular cluster.

How do I create a csv file containing the cluster label and its representative sequence?

(note: a cluster that appeared once in this file can re-appear somewhere in the middle of the file even though it is not shown in this case)

Thank you in advance for the help.

Answer: I created my own answer for this question.

```
# this program takes in cltr file and create a mapping file
#cluster,representative sequence name

import re
# Using readlines()
file1 = open('sub_SRR5892231_90.clstr', 'r')
file2 = open('map.sub_SRR5892231_90.clstr.csv', 'w')
Lines = file1.readlines()

# since I am using a while loop I need the following index
LineIndex = 0
while LineIndex < len(Lines):
# check if >Cluster lines are there
FF = re.search(r'^(>Cluster.*)', Lines[LineIndex])
# if the lines re found, then prep it for first column
if FF != None:
#print(FF.group())
firstcol = FF.group().replace('>','')
file2.write(firstcol+',')
print(firstcol, end =",")
LineIndex = LineIndex + 1
while LineIndex < len(Lines):
KK = re.search(r'^(.*\*)', Lines[LineIndex])
if KK != None:

seccol = KK.group()
seccol = re.sub(r'.*>','',seccol)
seccol = re.sub(r'\.\.\. \*','',seccol)
print(seccol)
file2.write(seccol+'\n')
break
LineIndex = LineIndex + 1
LineIndex = LineIndex + 1

file1.close()
file2.close()

```
Question: This might be a naive question, sorry.
How are "hub genes" determined in a PPI network? Just determine the nodes degree and choose top 10 or so?

Answer: A hub is a node that has a higher degree than other nodes in the graph. How this translates in practice seems to vary between papers. Some just use an arbitrary threshold either as a fixed number or as a percentage, others choose a threshold based on some statistical argument e.g. outliers in the distribution when assuming a scale-free graph.
As a side note, be careful about what you say about high-degree genes/proteins. Many of them have high degrees simply because they have been more studied (e.g. the historical oncogenes such as the RAS family, SRC, EGFR...) or because they belong to a particular class (e.g. ubiquitin-like families, chaperones) or because they are sticky proteins (some are well documented false-positives such as ferritin and ribosomal proteins for yeast two-hybrid or collagens in AP/MS, see also for example the [crapome database][1] for "contaminants" in AP/MS).
Finally, note that the [scale-free assumption probably doesn't hold for many biological networks][2].

[1]: http://reprint-apms.org/
[2]: https://arxiv.org/abs/1801.03400
Question: Hi,

Is there a way to do 3d plots in EIGENSOFT, using the ploteig.perl program? I have the 2d plot of the first two principal components from EIGENSTRAT, but I need to plot the first 3 principal components in a 3d plot. I was wondering if there are other programs/software which can do this?

Any help would be much appreciated.

Answer: You can use R with the scatterplot3d library to make a 3Dplot. Input file is the Eigensoft .evec file, which contains the first 10 principle components (I named this file vcf2eigen.pca.evec).

You can see a picture of the plot here: http://rpubs.com/rbagnall/89838

Here's the code:

```r
library(scatterplot3d)

# read in the vcf2eigen.pca.evec file
df <- read.table("vcf2eigen.pca.evec")

# give column names
names(df) <- c("SampleID", "PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10", "POPULATION")

# plot
with(df, scatterplot3d(PC1, PC2, PC3, color = as.numeric(POPULATION), pch=19, main="PCA 1000 Genomes Exome Data (3D)"))

# add legend
legend("topleft", pch=19, col=c("black", "green", "red", "blue"), legend=c("AFR", "EUR", "SAS", "EAS"))
```
Question: Hi All,

I am trying to make cnet plots of differentially expressed genes between 2 tissue types, looking at specific pathways identified from enrichGO analysis. When I plot this a bunch of non-significantly regulated genes are appearing as uncolored dots and I am not sure how to stop this from happening so that my plots only include significantly regulated genes within the pathway of interest. I have attached an example of one of the cnet plots generated and am hoping someone knows how to fix this problem.

Here is the code I am using for the cnet:

>selected_pathways = "response to oxidative stress"

>cnetplot(eGO2, foldChange = geneListC57mnx_chow_iWAT_vs_eWAT, showCategory = selected_pathways, circular = FALSE, colorEdge = TRUE)

And the code for generating eGO2:

>eGO2 <- enrichGO(gene = as.character(C57mnx_chow_sorted_gene_list$ensembl_gene_id),
OrgDb = org.Mm.eg.db,
keyType = "ENSEMBL",
ont = "ALL",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)

The file geneListC57mnx_chow_iWAT_vs_eWAT only includes genes that have a p.adjusted > 0.05 so I am confused as to why the plot is such a mess of insignificant genes. Any help would be much appreciated!

![enter image description here][1]

[1]: /media/images/bebf9529-1b09-4775-98a9-fb7449f9

Answer: I think the issue is that you're creating `eGO2` with non-significant genes. Try subsetting `C57mnx_chow_sorted_gene_list` to your thresholds before sending it through `enrichGO`
Question: The plot is from this publication, see page 2 the bottom of the plot with ACTG letters in different colours and sizes.

Could anyone clarify in layman's terms what it represents?

More importantly, how do we plot it using R or any other tools?

Here is the plot, if you can't access the paper:

Found related post at SO: Plotting a “sequence logo” using ggplot2?

Answer: The plot shows the comparison of the reference human sequence with N-number of samples (patients). The samples turn to have variation of sequences in that region. Now from these samples, you can compare the probability how often a given nucleotide is mismatched/mutated in sapmles. If the letter is big it means the sequence is invariant.

You can use weblogo either in their webpage or in R to make that figure.

Or if you have Chip-seq peaks, you can do motif analysis and make those plots using HOMER tools.
Question: I am learning how to analyse NGS data. I have a data for 192 samples. These were obtained through a targeted sequencing library prep.

I have 192 samples, but technically I have received 2 files for each sample. For example:

- sample1_TTGCCTT_L008_R1_001.fastq.gz
- sample1_TTGCCTT_L008_R2_001.fastq.gz

Presumably the reason there are 2 files is because paired-end sequencing was performed. I've been reading around on the various steps of NGS analysis but I can't seem to find an answer to the following question:

How to handle paired-end reads? Do you do have to merge them and if so when? Before alignment presumably? Also, do I have to uncompress the fastq.gz before I do anything? I am very new to NGS so apologies if these are really basic questions. Thanks.

Answer: You don't need to merge the R1/R2 reads. You don't say what kind of data this is but generally if you are aligning to a reference then you would use the two files together with an NGS aligner. Since the files contain reads from the same fragment their alignment to a reference provides spatial information.

All extant NGS tools should understand gzipped files. You should not need to decompress then during analysis (note: there may be some exceptions depending on very specific programs).

(An aside: If the reads are longer than the 1/2 size of the insert then they can overlap in the middle. )

Reads will overlap in this case

|------------------------------>100 bp| R1 - 150 bp
|-------------------------------------| Fragment 250 bp
|100 bp<------------------------------| R2 - 150 bp

and will not here

|--------> | R1 - 100 bp
|-------------------------------------| Fragment 350 bp
| <---------| R2 - 100 bp
Question: I know this question has been asked and answered many times, but I'm having trouble with making a gene function annotation on a series of differentially expressed genes for my tomato transcriptomes (Solanum Lycopersicum).

Our genomics group sent to me a series of top differentially expressed genes in my tomato experiment, with each gene having an annotation from the sol genomics network. The gene IDs look like this:

```
Solyc11g020960
Solyc04g081710
Solyc07g043000
```

I just need to classify them by their general function. I already have their protein products annotated in my dataset, I need to see what function and pathways they belong to. Naturally I started using DAVID but it doesn't seem to recognize these Gene IDs regardless of what identifier I choose. Also DAVID hasn't been updated in the last 5 years, which is also concerning. Tomato plants are a very well annotated organism, but I can't seem to find any online gene function annotation tools which have their metabolic pathways.

Any advice on how to get these gene functions annotated in a user-friendly way? The end goal is to find out what pathways these genes belong to and what functions they control, and display these findings in a useful way where I can read them.

Thank you and have a great day!

Answer: You can try PANTHER (http://www.pantherdb.org/). It can identify Solanum Lycopersicum.
Question: The Genbank ftp site (ftp://ftp.ncbi.nih.gov/genomes/) used to contain a folder called Bacteria, with all the bacterial genomes, but it seems to have disappeared. Does anyone know where these might have been moved to? I can't find anything about it on the NCBI site, and I need to access some gff files ASAP.

Answer: Assembled genome sequence and annotation data for GenBank genome assemblies is now (from 02-DEC-2015) available under: ftp://ftp.ncbi.nlm.nih.gov/genomes/genbank/
Question: Hi all, I have just a very simple question about subsetting data.

I have a data frame with many columns, two columns are indicating the fold change after RNA-seq data analysis.
This two columns contain positive and negative values. I would like to keep the opposite values (negative in one column and positive in the other and vice versa) while removing those which are the same sign (negative and negative in the two columns or positive and positive).

Example:

df$1 df$2
1. 2
-1. 1
-1. 4
3. -2
6. -4
-6. -6
6. 1
-3. 7
......

Resulting dataframe :

df$1. df$2

-1. 1
-1. 4
3. -2
6. -4
-3. 7
......

Any help, highly appreciated, thanks!

Answer: You can define a new column:

df$3 = df$1 * df$2

then

df <- subset(df, df$3 < 0, select = 1:2)
Question: Hello

I was looking through my created fasta files and realized that the number of bases are consistently 6 above the number that it should be according to the coordinates. I am not a strong scripter so how would I remove the last 6 bases from every sequence in my fasta file? Note: each line in the sequence is on a new line so I can't just remove 6 from the end of every other, I would need to iterate through it as a fasta file format.

>gene1
CACTGATTTGAGTTTTTTTCATAAATCAGAAACCGTTTGAATTATAAAAA
AAAAAACCTCACACAGCTCAAACTCAACCTCTGAAAAATAACCCCCAAAG
TGTAGCACTTCTAGCCCAATTCTTTAAACTTTGGTTGAAGGCTTCTGCAT
AGAGACGCGGGAAAGACAGTTTTACTGTTTAGCACTCTATGGAGCAACAT
CTGTAGCAACACTACTGGGGGGCCAGCGCGAGTATATGGACACAAAACAT
CATGTAGTGTAGGATTTCTGAATAGCAATACACCCTTTGTGGTGATGTAA
CAATAAAGGAAAGGGCATATTTTTGATGATCATGAGGTGTAGCCCCT

I would like to remove "GCCCCT" from the end. This is a multifasta file so there are hundreds of similar cases within the same fasta file. Thank you!

Answer: If you have the gene name, and they match the header in the fasta file, it is still fairly trivial to use python or biopython to trim your sequences.

1. Linearize your fasta file from multi-line to single line. Pierre's awk statement is the easiest option in the link posted in my comment above. This will make it easier to work with.
2. Make a file, tab-delimited with the gene name, start coordinate, and end coordinate
3. Save the code below as `trim_fasta.py`, and run as `trim_fasta.py coords.txt input.fasta output.fasta`

See code below:

#!/usr/bin/env python
import sys

with open(sys.argv[1], 'r') as f:
coords = {}
for line in f:
coords[line.strip().split('\t')[0]] = (int(line.strip().split('\t')[1]), int(line.strip().split('\t')[2]))

with open(sys.argv[2], 'r') as fasta:
with open(sys.argv[3], 'w') as out:
for line in fasta:
if line.startswith('>'):
header = line.strip().split('>')[1]
out.write('>' + header + '\n' + next(fasta).strip()[coords[header][0]:coords[header[1]])
Question: AFIK, Short sequencing and short pairwise sequencing generated by next-generation-sequencing is common now-a-days and there exists several alignment tools for them.

But I want to ask about long pairwise sequencing. There exists some alignment tools for long sequences like [rHAT][1]. That means there is at least one tool which generates long reads.

Now, I am asking that is there any tool which generates **Long Pairwise Reads**? if yes, then is there exist any aligner or mapper for that? Give me some idea about it.

Thanks in advance.

[1]: https://github.com/derekguan/rHAT

Answer: Depends on your definition of "long". The idea of paired end sequencing originally arose with clone based sanger sequencing when sequencing primers would be designed to read in from both ends of the insert. This was used with short inserts (such as a standard plasmid based cloning vector) to a) effectively extend the maximum read length, OR provide two readings of a sequence for error checking. With long insert libraries, such as BACs, paired-end sequencing allowed placement of sequence on scaffolds, and detection of structural rearrangements and copy number alterations.
Question: Hello all,

I have some RNA-seq data from an unusual biological system, and as a result I need to look for genes that are differentially expressed genes between two samples, one of which was *super* stressed. My DESeq2 results said that more than 70% of genes had an adjusted p-value below 0.1! I was expecting a lot of differentially expressed genes, but this seems very high. Is there such a thing as 'too many' DEGs? Is there a way I can test whether my data fall into this category? Does anyone know of any solutions to this issue? Thanks in advance!

Answer: This question has been asked before, see previous discussion at [Question: Differential expression for two very different samples][1]. As [swbarnes2][2] said, one assumption of DESeq2 and edgeR is that not too many genes are differentially expressed. edgeR allows the user to mitigate this issue with the argument `logratioTrim` in the function `calcNormFactors()` - see [here][3]. I don't know DESeq2 has a similar option - I believe it doesn't.

Regarding cut-offs: the author of DESeq2 argues that instead of applying a cut-off after a regular test for differential expression (that is, a test against the null hypothesis of zero LFC), it is better to test if the LFC is above a certain threshold. He argues that the later approach keeps the nominal fdr, while the former doens't. See discussions on DESeq2 vignette and [paper][4], and at the BioConductor forum [here][5] and here. Use the arguments `lfcThreshold` and `altHypothesis` from `results()` to set the desired LFC threshold and test.

[1]: https://www.biostars.org/p/262676/
[2]: https://www.biostars.org/u/4407/
[3]: https://support.bioconductor.org/p/94661/
[4]: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8
[5]: https://support.bioconductor.org/p/85509/
Question: Hi,

I downloaded copy-number-variation data from the TCGA portal. However, I found that for some samples could be present multiple replicates. Supposing I want a single copy number variation profile for each sample, how can I select the best one or aggregate the replicates?

Thanks

Answer: - Considering samples from the same patient 'replicates' is potentially dangerous. Spatial heterogeneity in solid tumors may mean that two samples have very different profiles.
- If the primary difference is just multiple platforms from the same tumor subsection, you can probably just use the highest-resolution method and call it a day.
- If you want to get complicated, you can look at this method out of Terry Speed's lab, but my feeling is that it's probably more hassle than it's worth, in general. http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2660872/
Question: Hi,

Would you please recommend me a good software for searching domains in a protein or coding nucleotide sequence? (except NCBI conserved domain search). Your help is greatly appreciated.

Answer: My favorites for proteins are

- [InterPro][1] (where you can find domains and get lots of information using InterProScan). You can also use InterProScan with nucleic acids
- The [UniProt database][2]

[1]: http://www.ebi.ac.uk/interpro/
[2]: http://www.uniprot.org/