cannin/biostars_qa huggingface dataset:
Question: Hi,

I am computing the lambda statistics in R to check whether the distribution of my p-values follow a null distribution using:

pvalue <- runif(1000, min=0, max=1)
chisq <- qchisq(1 - pvalue, 1)
lambda <- median(chisq) / qchisq(0.5, 1)

Is there a procedure for constructing the confidence interval for this lambda genomic factor in R?

Thanks

Answer: I encourage you to follow up on my points here with a statistician in your department. Always better to get more brains looking at the same thing. There are various ways to calculate confidence intervals for series of data.

You could start with a typical formula for the standard error of the mean (SE mean), which is the standard deviation divided by the square root of the sample n:

σ / sqrt (n)

Lambda is a median value, so, we have to convert the SE mean to SE median:

1.253 * (σ / sqrt (n))

(various articles on the WWW validate this)

As you are dealing with a Chi-squared quantile distribution, in order to calculate the confidence intervals, I believe you can do this with the following:

pvalue <- runif(1000000, min=0, max=1)
chisq <- qchisq(1 - pvalue, 1)
lambda <- median(chisq) / qchisq(0.5, 1)

lambda
[1] 0.9963626

SE.median <- qchisq(0.975, 1) * (1.253 * ( sd(chisq) / sqrt( length(chisq) ) ) )

SE.median
[1] 0.008914307

lambda + (SE.median / qchisq(0.5, 1))
[1] 1.015957

lambda - (SE.median / qchisq(0.5, 1))
[1] 0.976768

I increased the number of random p-values to 1 million just to make the confidence intervals tighter. At 1000 (your original post), they were pretty wide.

Take this as a first step into this. I encourage you to look up all functions / formulae that I have used.

Kevin
Question: I would like to do RNA-seq differential expression analysis on an experiment run on 8 pigs (4 cases and 4 controls, not matched). I've already done the tophat/cuffdiff approach using the pig genome. My PI also wants me to try out Trinity (without using the pig genome sequence).

I noticed that trinity can take in multiple fastq files at a time. How should I run my Trinity run(s)? Do I do 8 separate Trinity runs? Two runs for the two conditions?

Thanks

Answer: If you're looking to do a DE analysis on your Trinity results, then provide all of your fastq's at the same time. Take note if you have PE reads or SE reads. They are supplied differently to Trinity on the command line. This will create a single *de novo* assembly. If you assemble each sample separately, you'll have trouble with the RSEM step when aligning the raw reads back to the assembly to get abundance estimates. You can do individual runs, if you simply just want a *de novo* assembly for that sample. But for the DE analysis, you'll need to create a single assembly, combining all samples as inputs.

Example:

if PE: `cat` all your read 1's, and `cat` all your read 2's

Trinity --seqType fq --JM 50G --left reads_1.fq --right reads_2.fq --CPU 6

For the DE analysis downstream you'll be able to specify replicates, and sample names.

If all fastq's are SE there is a `--single` option to use instead of specifying read orientation.

Check out [the Trinity site](http://trinityrnaseq.sourceforge.net/#typical_options); and for full usage there is a command line parameter: `--show_full_usage`

Hope this helps.
Question: Hi Biostars Community,

I have searched the forum, but couldn't find a perfect answer for this.

I have downloaded some RNA-seq data from GEO ([GSE112656][1]). They are basically 11 osteoarthritis samples and 11 rheumatoid arthritis (so 22 in total).

I also have 4 samples of RNA-seq data generated from my lab with three technical replicates (so 12 total).

1. Drug A - treated
2. Drug A - control (untreated)
3. Drug B - treated
4. Drug B - control (untreated)

I want to use the GEO ([GSE112656][1]) and also my lab data to conduct PCA to see which ones are similar and cluster together and which ones are different.

How would I normalize and transform these datasets before following these directions:

https://bioconductor.org/packages/devel/bioc/vignettes/PCAtools/inst/doc/PCAtools.html

I was considering following this post:
https://www.biostars.org/p/411830/

Basically, I am planning to combine all the samples (22 +12 = 32 samples) into one large data frame and then generate log2 transformed TMM followed by throwing that table into a PCA function to visualize the PC_1 by PC_2 table, or maybe PC_1 by PC_2 by PC_3 to seeing how they cluster/how similar or different they are to one another. I am planning TMM because [HBC][2] recommends this to compare between samples and within samples. Would all this be best practices?

Please, if you can suggest something better, I would appreciate it.

[1]: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE112656
[2]: https://hbctraining.github.io/Training-modules/planning_successful_rnaseq/lessons/sample_level_QC.html

Answer: Hey mate,

I would process both datasets independently to produce variance-stabilised expression levels via *DESeq2*, bind these together, and then eliminate the [assumed] batch effect via `limma::removeBatchEffect()`. It is the resulting batch-corrected, column-bound, variance-stabilised expression matrix on which I would conduct PCA, with `scale = FALSE` (but `center = TRUE`).

No issue going the *EdgeR* / *limma*-*voom* routes, either, in which case use `log2 CPM+1` in place of variance-stabilised expression levels.

Kevin

# @@
Thanks for up-votes and accepting the answer. It would also be interesting to process all samples together from the beginning, and then still remove `batch` on the variance-stabilised expression matrix that is produced - see https://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#why-after-vst-are-there-still-batches-in-the-pca-plot
Question: Hi,
I read the GSE file into R as follows.

gse <- getGEO("GSE4928", GSEMatrix=TRUE)
Now I want to convert all probe IDs to gene symbols and write the entire annotated dataset back to my computer. What should I do?
By the way, which dataset is stored in gse? SOFT formatted family file, MINiML formatted family file or Series Matrix File?

Answer: Take a look at this code:

gse = gse[[1]] # get just the first element in the list
head(fData(gse))
symbols = fData(gse)[,'Gene Symbol']

Now, you have the gene symbols. The gse object (after choosing the first element in the list using gse[[1]]) contains the information from the Series Matrix file AND the GPL file.
Question: Hello,

I have a FASTA file, and I want to extract a subsequence from it.

```
>first
AGCATGCATCGTACGATCGTACGATCATCATC
GTACGATCGATCGATCGATCGATCGTGTGCTA
GACGATCGATCGATCGA
>second
GTACGATGCATCGATCGATCGATCGATCGATC
GTGATCGTGTCGTCGTCGTCGATAAGCATCGT
AGTACGATCGATCGTA
```

I know I can do this with samtools:

samtools faidx seq.fa first:1-5

output is:

```
>first:1-5
AGCAT
```

However, my actual sequence names are much more complicated, and anything other than one word doesn't seem to work.

If `seq2.fa` is:

```
>first seq
AGCATGCATCGTACGATCGTACGATCATCATC
GTACGATCGATCGATCGATCGATCGTGTGCTA
GACGATCGATCGATCGA
>second seq
GTACGATGCATCGATCGATCGATCGATCGATC
GTGATCGTGTCGTCGTCGTCGATAAGCATCGT
AGTACGATCGATCGTA
```

Then the output of:

```
samtools faidx seq2.fa first:1-5

>first:1-5
seqAG
```

​```​
samtools faidx seq2.fa first seq:1-5

>first
seqAGCATGCATCGTACGATCGTACGATCATCATCGTACGATCGATCGATCGATCGATCG
TGTGCTAGACGATCGATCGAT
>seq:1-5
```

```
samtools faidx seq2.fa 'first seq':1-5

>first seq:1-5
```

Also, if my sequence description contains a number (like `seq3.fa`):

```
>seq123
AGCATGCATCGTACGATCGTACGATCATCATC
GTACGATCGATCGATCGATCGATCGTGTGCTA
GACGATCGATCGATCGA
>seq456
GTACGATGCATCGATCGATCGATCGATCGATC
GTGATCGTGTCGTCGTCGTCGATAAGCATCGT
AGTACGATCGATCGTA
```

```
samtools faidx seq3.fa seq123:1-5

>seq123:1-5
```

Is there a way to use samtools to extract a subsequence when the descriptor is more complicated than just a word? I've tried with samtools 0.1.18 and 0.1.19.

Thanks

Answer: faidx stops reading the name after the 1st blank: https://github.com/lh3/samtools/blob/master/faidx.c#L89

while ((ret = razf_read(rz, &c, 1)) != 0 && !isspace(c)) {

you could fix this by converting the spaces to an underline

tr " " "_" < in.fa > out.fa && samtools faidx out.fa
Question: As part of their submission requirements ENA/ArrayExpress require BAMs to include all reads (e.g. also unmapped ones). TopHat in its wisdom splits unmapped reads from the accepted_hits.bam output.

However, creating a 'complete' BAM is not so straightforward due to problems in the unmapped reads BAM (detailed on seqanswers). I've tried the script mention on SeqAnswers, but it doesn't solve all the problems for me. I still get unmapped reads which don't have a mate in accepted_hits.bam

Given TopHat's popularity this should be a solved problem, right? So, how does everyone deposit their TopHat alignments for publication? Or are fastqs the preferred format.

Any suggestions?

Answer: With the help of Christian Brueffer, I managed to resolve this problem. See his github: https://github.com/cbrueffer/tophat-recondition/issues/1#issuecomment-129849801
Question: Hi All,

Do we have timeline for human hg39 genome? I ask these questions because I have lots of hg37 data need to transfer to hg38, but I am wondering if hg39 come out soon, then I think I can wait for hg39 version.

Thanks

Answer: https://www.ncbi.nlm.nih.gov/grc/human

> the GRC remains committed to its mission to improve the human reference genome assembly, correcting errors and adding sequence to ensure it provides the best representation of the human genome to meet basic and clinical research needs. We will continue to make these updates publicly available at regular intervals in the form of patch releases, but **have decided to indefinitely postpone** our next coordinate-changing update (**GRCh39**) while we evaluate new models and sequence content for the human reference assembly currently in development.
Question: Hello everyone,
I have a huge dataset with a bunch of human samples to analyse. Of course, I run into troubles because the samples come from different donors and when I PCA those samples, well...it's a bit dodgy. They cluster according to their condition but I am not sure about how am I supposed to deal with this batch effect? A few time ago, I used SVA package but I wasn't happy with that.

A problem related to this is probably due to the fact that my samples are not trimmed appropriately. I have a lot of problem with the facility that generated these fastq files because sometimes they provide me trimmed samples, sometimes they don't (given the fact that this whole dataset comes from different batches/years).
Thus, my questions:

1. Don't you think that all of my samples, to generate useful data, must have been processed in the same identical way (e.g. same Sliding window, leading, trailing, minlen)? I am quite confused about this.
2. What if, by any chance, I trim an already-trimmed file?
3. When I am trying to trim my samples, I don't manage to remove adapter contamination..according to my beloved multiqc report there's a huge nextera transposase sequence contamination that Trimmomatic can't remove, even when selecting specific adapters...

Yours,
M

Answer: > A problem related to this is probably due to the fact that my samples
> are not trimmed appropriately.

I wouldn't be so sure. For instance, I know that STAR aligner is pretty robust to having wrong sequence on the ends of reads.

If your trimmer isn't trimming anything, maybe nothing needs trimming. If you have a big batch effect, that's likely real, and not an artifact you can fix.
Question: In looking for a way to check if a bamfile is truncated I noticed that this does the job very quickly.

samtools view file.bam | 1>/dev/null

The thing is I typed the pipe by accident; the command doesn't make sense to me but when you remove the pipe the command takes much longer (Im not sure it ever finishes). My question is why does this work/what is going on.

Answer: When you run without the pipe "|" symbol

samtools view file.bam 1>/dev/null

1 is interpreted as stdout so it is exactly the command we all are very used to:

samtools view file.bam > /dev/null

You know it. This is a way to read bam file without header into a SAM file if you replace "/dev/null" with "your.SAM". It takes very long time. /dev/null just discards all the data that was put into it, put it actualy first has to get it from the samtools. The whole file is read despite this /dev/null anyway.

When you add the "|" pipe symbol the way you mentioned:

samtools view file.bam | 1>/dev/null

Your first command tries to read the whole bam file without the header and pipe it (kind of sam file) via stdout into a second command, in your case it is "1". Most likely you do not have program or script called "1" in your system. So the second command will end immediately with an error and as it is the last command in your pipe before redirect it will terminate the whole pipe. What you see on your screen is what the first command had time to write as an error message into stderr since stdout was redirected with pipe. When bam file is truncated stderr contains information that file was truncated when it is a valid bam file stderr is empty so no message. You can test all this by running this:

samtools view file.bam | biostars

since you most likely do not have a program called "biostars" in your system it will behave almost as your "samtools view file.bam | 1>/dev/null" other than it will clearly tell you that there was no "biostars" command found. Then you can run command

samtools view file.bam 2 > err.txt

terminate it immediately with Control+C since it will try to print the whole sam file into your stdout. Now take a look at the contents of the err.txt. It will be empty for proper bam and will have an error message for truncated one. This means that samtools already raised an error before you terminated it manually a moment after it started to read a bam file. So if you do

samtools view file.bam 2 > err.txt | wow

you will get information whether bam is normal or truncated in err.txt file if you do not have program called "wow" in your system =)
Question: I am looking for a **clear** explanation on what WGCNA actually produces as output.
I've seen many papers using it and **I actually read both the manual and the paper**, but I do still not really understand what exactly these hubs and hubgenes and modules are and in which situation one would use WGCNA.

Could someone please explain for dummies what results this tool actually returns and most importantly, if this is possible, an example in which situation a researcher would use it, so what the biological question could be that motivates to use use WGCNA. Thanks!

Answer: The idea behind WGCNA is simply to view correlations between gene expression profiles as edge weigths of a graph. That is the correlation matrix is viewed as the adjacency matrix of a graph in which the nodes are genes and the edge weights are the correlations. From there one could then apply various graph algorithms to characterize the graph. What the output is depends on which algorithm is applied. Modules are another name for clusters and are obtained by applying a graph-based clustering algorithm. Hubs are nodes that have a high degree (i.e. are highly connected).
As to what questions can motivate the use of WGCNA, there are many. For example you could be interested in finding groups of genes that are co-regulated. These might be detected as modules in the graph. If interested in identifying "master regulators" of a process, you may find them as hubs. For examples of use, I suggest to check the literature. [Here][1] is one random example and [here][2] another one.
Note that the methods in WGCNA are generic and any similarity matrix could in principle be used as input, not just correlations and also that there are other graph algorithms than those available in the WGCNA package.

[1]: https://www.frontiersin.org/articles/10.3389/fgene.2018.00508/full
[2]: https://www.nature.com/articles/s41598-020-66132-4
Question: Does anyone know if there is a glycomics database with info on sugar modifications etc of proteins?

I have a few protein/gene names and would just like to get some indication of whether they're known to be glycated and whether the sugars have been identified.

Answer: Answering my own question (but not in a remotely exhaustive or even particularly thorough manner 😅):

- https://unicarb-db.expasy.org/
- https://acgg.asia/db/gpdb/
Question: Dear All,

I have two sets of annotation files for 15 Bacterial genomes. One set from NCBI annotations (from RefSeq) and the other from Prokka (I have run it in my local machine).

Which one is advisable to use between the two for all downstream analysis?

Awaiting your valuable feedback

Thank You

Answer: If I remember correctly, prokka comes only with HAMAP database of HMMs, which will produce terrible annotations on prokaryotic genomes. To get good annotations you would need to install at least Pfam and TIGRfams. Don't know if you have done that or not, but you can find out by looking at prokka's annotations. If there are many hypothetical proteins for prokka where NCBI files have meaningful annotations, chances are that you don't have any extra prokka HMM databases. If you are literally comparing identical genomes, it may be better to go with NCBI annotations.
Question: Dear list,

I'm trying to cluster data and then plotting a heatmap using heatmap.2 from ggplot, the script works perfectly with matrices up to 30000 rows, the problems is that I'm using matrices up to 500000 rows (data_sel), and when I try to cluster I get this error:

heatmap.2(as.matrix(data_sel),col=greenred(10), trace="none",cexRow=0.3, cexCol=0.3, ColSideColors=fenot.colour, margins=c(20,1), labCol="", labRow="",distfun=function(x) dist(x,method="manhattan"))
Error in vector("double", length) : vector size specified is too large

Is there any approximation using R to plot heatmpas with his big data?

Thanks in advance

Answer: With really large data, even a big-ass server may not be enough. The best advice I can offer is to try reducing the size of your data set. Doing this will require intelligently thinking about what it is you're trying to represent. (related: are you really going to be able to pick out 500k individual points on a little graph the size of your computer screen?)

If, for example, you're plotting expression data from an exon array, maybe you could merge the data and plot per-gene instead of per-exon. If you're looking for patterns of differential expression, maybe you could plot just the top 1000 most differentially expressed. (and so on)
Question: As part of a large-scale genomic integration project, I ended up with chromosome positions mapped to human genome patches like:

- `HG104_HG975_PATCH`
- `HG1257_PATCH`
- `HG19_PATCH`
- `HSCHR10_1_CTG2`
- `HSCHR17_1`

I need to map these patch IDs back to chromosome numbers to unify them with rest of the database. While the last two cases indicates chromosome numbers as 10 and 17 respectively, am looking for a strategy to map the first three.

Full list of patch IDs that we have is uploaded [here][1].

Thoughts?

[1]: http://pastebin.com/zGmWPmsd

Answer: Here you go.

Note that Ensembl version 75 is based on the GRCh37.14 assembly. Version 76, expected later this month, will be based on GRCh38.

```
mysql -h ensembldb.ensembl.org -u anonymous
Welcome to the MySQL monitor. Commands end with ; or \g.
Your MySQL connection id is 659408
Server version: 5.1.72-log MySQL Community Server (GPL)

Copyright (c) 2000, 2013, Oracle and/or its affiliates. All rights reserved.

Oracle is a registered trademark of Oracle Corporation and/or its
affiliates. Other names may be trademarks of their respective
owners.

Type 'help;' or '\h' for help. Type '\c' to clear the current input statement.

mysql> use homo_sapiens_core_75_37;
Reading table information for completion of table and column names
You can turn off this feature to get a quicker startup with -A

Database changed
mysql> SELECT DISTINCT s1.name AS patch, s2.name AS chromosome, a.exc_seq_region_start AS start, a.exc_seq_region_end AS end FROM assembly_exception AS a, seq_region AS s1, seq_region AS s2 WHERE a.seq_region_id = s1.seq_region_id AND a.exc_seq_region_id = s2.seq_region_id AND a.exc_type IN ('PATCH_FIX', 'PATCH_NOVEL');
```

https://gist.github.com/RamRS/82dfc8414c6c83853c544226316eeed2
Question: Hi,

In various CNV analysis papers there is mention of excluding areas of the Genome from analysis

e.g. if the CNV spanned a centromere, if the CNVs overlapped at least 50% with regions previously defined as being prone to false positives due to somatic mutations, etc

Would anyone have a paper that details a list of which regions to exclude from CNV analysis, or a preprepared file for PLINK?

Best regards

Bob

Answer: There are several tracks in UCSC that might be useful for selecting regions to exclude from CNV analysis. For example if you wanted to exclude near centromeres and telomeres, regions of high or low GC, assembly gaps, regions with poor mappability, repeat regions, etc. Here are some examples of such tracks with links to their descriptions

Chromosome Band

Gap Locations

GC Percent

GRC Incident

HiSeq Depth

Mappability

Repeat Masker
Question: How to extract features from protein sequences, so that it can be converted into vector for training the data in machine learning. From some papers I found methods like using AAindex, PSSM for training data. But I was unable to find the detailed method behind it. Please, suggest some papers or links which can be helpful.

Answer: Features you want to extract are divided into two groups; 1) sequence based sequence, 2) features extracted from the predicted structure.

Amino acid composition, amino acid property, amino acid distribution and etc are in group one. There are mainly two R packages, *Seqinr* and *BioSeqClass* from Bioconductor. I attach [a table from my thesis][1], which summarize this threat.

I would recommend you reading this paper; and the code is implemented in the BioSeqClass package.

Prediction of protein folding class using global description of amino acid sequence. PNAS, 92(19):8700-8704, 1995 (BioSeqClass - Bioconductor package)

[1]: https://drive.google.com/file/d/0B1CZNSfA7sV8X0hyYlFKTzRZbUU/edit?usp=sharing
Question: Hi,

I am generating overlaps using this command in bedtools:

intersectBed -a snps.bed -b region.bed -wb

And now I am trying to compare the same files using the R package GenomicRanges, countOverlaps(snps, region.bed, type=c("within"), but I am not getting the same results.

In particular, it seems as the positions at the borderline of the region.bed are not counted as overlaps,
for example this entry in snps.bed:
1	112162878	112162879

and this entry in region.bed:
1 112162834 112162878

Is this what others have seen as well? How can I get bedtools to count these entries as overlaps, without adding flanking bases?
Thank you

Answer: BED files are generally 0-indexed, so you're generally not going to get overlaps with applications that do operations on input files that use 0-based indexing, which generally do not count endpoints in calling overlaps. GRanges documentation says it uses 1-based indexing, so it will include endpoints in generating overlaps.

Your solution is perhaps not to count those kinds of overlaps in BED file inputs using GRanges, by specifying that your (BED-formatted) inputs are 0-based and need adjustment in R. The [documentation][1] for GRanges explains this in more detail; see the `starts.in.df.are.0based` option.

[1]: https://bioconductor.org/packages/3.3/bioc/manuals/GenomicRanges/man/GenomicRanges.pdf
Question: I am trying to include the --keep and --remove options in my plink2 command. I am finding that despite my files for these options having identical text to the main file's IDs, all the samples are removed.

Command: plink2 --bfile pca_final --keep EUR.sample --remove toberemoved.txt \
--make-bed --out unrelated_EUR

So my original .fam file has:

0002

The --keep file would have:

#FID	IID
0002	0002

The --remove file would have:

#FID	IID
1067	1067

Answer: I was able to figure out the issue. The FID column (the first column in the fam file) was all zeroes and I put text there as I was confused by the format used by my colleague.