cannin/biostars_qa huggingface dataset:
Question: Hello,

I am analysing RNA-seq data to investigate differential gene expression in hybrids compared to parental species. Since I work with natural populations, I have few samples (5 of two different tissues for each species and hybrid).

I am using the DESEQ2 package for my expression analysis. What I observe is that many genes, although they have a high log fold change (more than 1 or even 1.5), they have a padj > 0.1. While this is true in one group, in another group, genes with a log fold change of 1 or even lower are having a padj < 0.1 or even padj < 0.05.

I was wondering what are the reasons for this observation?

Thank you in advance.

Answer: Fold changes tend to be higher when genes have overall lower expression (which means low counts).
Since low counts have lower power than high counts the significance for these fold changes is often low unless these FCs are supported by many replicates.

Example 1: Two genes had expression of 50 an 5. That would be a fold change of 10.

Example 2: Two genes had expression of 5000 an 500. That would be a fold change of 10 as well.

Still, the second one is much more reliable as the first one could be a product of the technical noise produced by the sequencing. Adding or reducing e.g. 10 counts to example a can change the result quite much:

`50 - 10` vs `5 + 5` would already change the original FC from 10 to 2.6 whereas

`5000 - 10` - `500 + 5` changes the FC from 10 to 9.88.

You can see that higher counts are less affected by small fluctuations in counts, therefore they are more reliable.
In DESeq2 you can check the `baseMean` column to get the average expression. This is probably low for many of these genes with high FCs but large padj. You can visualize this relationship of baseMean to logFC with the `plotMA` function.

This is where the concept of shrinkage kicks in. It aims to estimate the "true" fold changes from the data. As you can see below there is little evidence for the fold changes of the genes with low baseMean to be actually true, so they are shrunken towards zero. If you want lowly-expressed genes to be significant 8given they in fact are DEG) then you need most importantly many replicates and high sequencing depth.

Check the DESeq2 vignette for it.

Some examples:

Unshrunken FCs:

![enter image description here][1]

Shrunken:

![enter image description here][2]

[1]: https://i.imgur.com/8nsORxP.png
[2]: https://i.imgur.com/ufNK0sH.png
Question: Does anyone know of a simple scripting solution to take a list of accessions and pull out the gene name from the header of FASTA sequences?

For instance given the accession: XP_016469325.1

and given the FASTA entry:
>XP_016469325.1 Nicotiana tabacum|C3H|C3H family protein
MEEELLKRNTDCVYFLASPLTCKKGIECEYRHSEIARLNPRDCWYWLAESCLNPTCAFRH
PPLESHAETSSESAPPQHKSAVPVNKTNVPCYFYFNGYCIKGERCSFLHGPDDGTTTWKS
SKIASGVPDGPTAEKKTSVGSETGPASVEKPSNSSETGSKAAAHEYIKSQVDLISMTNDV
GEQSASHETSGSPSEEATAVRLDSLVPAEGFTQGGSDLSPDWSSDEEVEDNVEREEWLES
SPGFDVLVDDRIEGWSHKDDHSYLLQHDRECDERFAGYDFENNLEYDPAYPDMRIVSDEE
LDDSYYSKVESHEVNEYAREIVIPAHGRQSIPHKRKFPREPGFCARGNVDLRDLLKKRRV
IESDPPNYLSRRLDLSRFNAREQCRDRHRPQGSRWMPQSLASKLESNSSFSSGFVDATRL
EGANQLKKLRQSHRSSYRQQHFKDRRRGRSQPFANETPRRMASRQRSTEVPKIFGGPKTL
AQIREEKIKGREDGNSFERTVPSGGSEREDFSGPKPLSEILKDKRRLSSVVNFSN

I would like to have output the gene name "C3H" that is spanned by the "|"

This script I modified from a previous post can grab the gene names, however I'm not sure how to only get the gene names corresponding to a separate list of accessions (accessions.list).

with open('PlantTFDB_ALL_TF_pep.fas','r') as f:
for line in f:
if '>' in line:
line = line.strip().split('|')
print(line[1])

Answer: File `acc.txt` contains the accession IDs:

$ cat acc.txt
XP_016469325.1

Firstly, getting the mapping relationship between acc and gene name and saving them to acc2gene.tsv

$ grep '>' seqs.fa | \
perl -ne 'next unless />(\S+).+\|(.+?)\|/; print "$1\t$2\n";'
XP_016469325.1 C3H

Secondly, joining acc.txt and acc2gene.tsv using csvtk (http://bioinf.shenwei.me/csvtk/download/) or other tools：

$ csvtk join -H -t -k acc.txt acc2gene.tsv | csvtk cut -H -t -f 2

Original answer

Searching header line using `grep -f `, and capturing accession ID and gene name using regular expression in Perl.

$ grep '>' seqs.fa | \
grep -f acc.txt | \
perl -ne 'next unless />(\S+).+\|(.+?)\|/; print "$1\t$2\n";'
XP_016469325.1 C3H

If you just want the gene names, that would be easier:

$ grep '>' seqs.fa | grep -f acc.txt | cut -d '|' -f 2
C3H
Question: Hi there!
I´m using KIRC database from TCGA and downloaded FPKM values, then converted to TPM with this formulas:

FPKMtoTPM <- function(x) {
return(exp(log(x) - log(sum(x)) + log(1e6)))
}

df <- data %>%
mutate_if(is.numeric, FPKMtoTPM)

Now i want to convert TPM to log2(TPM+1). How can i do this?

Thanks!
David

Answer: Just as you defined the first function FPKMtoTPM you can define another one for the log transformation of the TPM values

logTPM <- function(x) {return(log2(x+1))}

df %>% mutate_if(is.numeric, logTPM)
Question: Hi All,

I tried to download some BS-Seq data from [McGill Epigenomics Mapping Centre (\[Here Database Link\]\[1\])][1]

I have a little assumption about the BS-Seq wig/bigwig files in http://epigenomesportal.ca/edcc/ database.

KNIH wig file is 0-10 scale, DEEP/CEEHRC wig files are 0-100 scale while blueprint and Roadmap are 0-1 scale. Anyone met this problem before?

It looks these wig file are coverage not methylation signal. the scale actually represent the sequence depth.

Thanks.

[1]: http://edcc-dev2.udes.genap.ca/cgi-bin/directory.cgi?b=2017-10&as=1&i=3&ctc=1&session

Answer: **I obtained a response from IHEC and it seems my assumption is right, each dataset is uploaded by different group.**

The annotation tracks deposited on the IHEC Data Portal have been generated by each respective data producing consortia, each with their own softwares and parameters. This explains the difference you see across tracks of different groups. IHEC is currently working on generating a pan-IHEC dataset using a unified data analysis pipeline and standardized quality control metrics. It will however take several more months before we reach this step. If you want to do real comparison across IHEC members datasets, I would currently recommend obtaining the raw data.
Question: Hello. I've been trying to use awk to solve the following problem for some time with no luck. Was hoping someone here could give me some clues as to what I'm doing wrong.

I have a file that looks like this:

INPUT:

>chr8:76290516-76290880
578	T
579	G
580	A
>chr14:22131464-22132025
468	T
469	G
470	A
>chr12:33695439-33695441
468	T
469	G
470	A

Each record in the file has a header that starts with `>` I would like to print a new column which is essentially a line number, starting after the header which starts with `>`, and to begin counting from the number following the `:` in the header. I would like each record (after the header `>`) in the file to be treated/counted independently.

I have tried to do this in steps, starting with adding the counts to the third column first, and then will attempt to add the header value to the counts following that. I have had no luck with getting the counts in the third column using this command below:

awk '{FS = "/n"}{RS = ">"}{if(!/^>/){print $1, $2, NF, $3 }}' input.txt

DESIRED OUTPUT:

>chr8:76290516-76290518
578	T	76290516
579	G	76290517
580	A	76290518
>chr14:22131464-22131466
468	T	22131464
469	G	22131465
470	A	22131466
>chr12:33695439-33695441
321	T	33695439
322	G	33695440
333	A	33695441

Any pointers would be very much appreciated. Thank you so much!!

EDIT: I did not include my failed attempt to solve this problem, as @mensur pointed out to me, so I have edited to include. Thank you.

Answer: awk -F '[:]' '/^>/ {print;P=$2;X=0;next;} {print $0, (P+X);X++;}' input.txt
Question: Dear all,

I've got a large list of variant as VCF format per chromosome (human), which annotated using VEP (release 98) as offline mode. But, I didn't get the corresponding rsid at the output. Could you please help me out how I can obtain the "rsid" during the annotation?

Thanks

Answer: [--check_existing=1](http://www.ensembl.org/info/docs/tools/vep/script/vep_options.html#existing)

> Checks for the existence of known variants that are co-located with
> your input. By default the alleles are compared and variants on an
> allele-specific basis - to compare only coordinates, use
> --no_check_alleles.
>
> Some databases may contain variants with unknown (null) alleles and
> these are included by default; to exclude them use
> --exclude_null_alleles.
>
> See this page for more details.
>
> Not used by default
Question: Hi everyone, I'm performing analysis of some RNAseq samples, and currently trying to cope with batch effect.plotting PCA of vsd transformed data, I can clearly see two batches which are differ fromt the others.

plotPCA(vsd, intgroup = c('batch'),ntop= 34085)
![PCA][1]

If I'm not mistaken, then for DE expression analysis I could use design formula of dds function, to reduce batch effect:

dds <- DESeqDataSetFromMatrix(countData = counts, colData = coldata, design =~ batch + type)
dds <- DESeq(dds)
vsd <- varianceStabilizingTransformation(dds, blind = TRUE)

But I also would like to visualise results of batch correction. For PCA plot I could simply follow the DESeq2 manual and just perform limma batch effect removal:

mat <- assay(vsd)
mat <- limma::removeBatchEffect(mat, vsd$batch)
assay(vsd) <- mat
plotPCA(vsd, intgroup = c('batch'),ntop= 34085)

Which gives me a nice results:
![batch corrected PCA][2]

But now I'd like to visualise expression of individual genes of interest. I could use vsd transforemd data, but these numbers mask the actual count numbers, therefore laking info on actual level of gene expression, which I'd like to preserve.
![expression of gene X fo given conditions, vsd transformed data, limma correction][3]

So what would be the better solution here? Could I use limma removeBatchEffect() on normalised count table, maybe?

Thanks!

[1]: https://i.ibb.co/bzspVz1/Intitial-pca.png
[2]: https://i.ibb.co/RQhzx2x/batch-corrected-PCA.png
[3]: https://i.ibb.co/Cs6jD0j/geneX.png

Answer: `vst` transforms to approximately log2-scale, that is simply a data transformation to allow a wide range of counts to fit on a narrow scale rather than having counts spread between 0 and numbers in the hundreds-of- thousands.
Question: I am attempting to learn a bit more about processing RNA-seq data, and I am starting off by trying to recreate some figures from a paper. To do this I need to import the data from GEO. They have stored they data in SubSeries and SuperSeries and I have been attempting to download the SubSeries and use the data. I initially tried importing that data into R using:

mRNA_data <- getGEO('GSE42379', GSEMatrix=TRUE)

When looking at the data, it downloaded into two separate GSE series. I have been unable to access the data from both of those GSE files, as it only seems to parse the first of the files. Is this because the data was obtained on two different platforms? How would I be able to access the data in the second GSE series? The only way I have found that will access all the data is if you use `GSEMatrix=FALSE` which parses all 28 samples but then to access them is a complete nightmare, which again is not all that well explained (or at least to my limited abilities/understanding). The data I am trying to access is [here][1].

Thanks

[1]: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42379

Answer: GEOquery always returns a `list` when using `GSEMatrix = TRUE` for exactly situations like this where the Series record has data from multiple platforms.

```
library(GEOquery)
mRNA_data <- getGEO('GSE42379', GSEMatrix=TRUE)
names(miRNA_data)
length(miRNA_data)
# mRNA_data is a `list`
human_data = mRNA_data[["GPL8759"]]
mouse_data = mRNA_data[["GPL8321"]]
human_data
```

Now, `human_data` and `mouse_data` are each `ExpressionSet` objects.
Question: Hi Senior Methylators:

This figure is derived from fastqc analysis to a BS-seq fastq data.

However, Why the content of "G" is lowest rather than "C"?

OK. Maybe I get the idea: it is from Sample_2.fastq?

Answer: Yup, you nailed it, it's because you're looking at sample #2, where you'll have a lack of Gs and an abundance of As.
Question: I'm trying to carry out an eQTL analysis on some TCGA BRCA data I have.

For the genotypes, I have a list of birdseed files, that each have a hybridisation reference that looks like: `WRIED_p_TCGA_143_147_150_Hahn_N_GenomeWideSNP_6_A09_799960`.

My RNA-Seq data is labelled with a portion of the sample barcodes, e.g., **TCGA.3C.AAAU** (which would correspond to the barcode **TCGA-3C-AAAU-01A-11R-A41B-07**).

For eQTL analysis, I need to pair up a genotype with the corresponding expression values, but I need a file providing the mappings. As per [the TCGA Website][1], I feel like its probably a SDRF file I need, but I'm not entirely sure about this, and I can't find one through the TCGA data portal in any case.

Could someone lend me some guidance?

[1]: https://wiki.nci.nih.gov/display/TCGA/Sample+and+Data+Relationship+Format

Answer: See if either of these help.

sdrf from Broad: https://tcga-data.nci.nih.gov/tcgafiles/ftp_auth/distro_ftpusers/anonymous/tumor/brca/cgcc/broad.mit.edu/genome_wide_snp_6/snp/broad.mit.edu_BRCA.Genome_Wide_SNP_6.mage-tab.1.2023.0/broad.mit.edu_BRCA.Genome_Wide_SNP_6.sdrf.txt

sdrf from UNC: https://tcga-data.nci.nih.gov/tcgafiles/ftp_auth/distro_ftpusers/anonymous/tumor/brca/cgcc/unc.edu/illuminahiseq_totalrnaseqv2/totalrnaseqv2/unc.edu_BRCA.IlluminaHiSeq_TotalRNASeqV2.mage-tab.1.1.0/unc.edu_BRCA.IlluminaHiSeq_TotalRNASeqV2.1.1.0.sdrf.txt
Question: In Seurat's pbmc3k tutorial, they set the CreateSeuratObject with various parameters including `min.cell=3` and `min.features=200`

Can anybody explain exactly what these two parameters mean and where they pull the numbers 3 and 200?

Answer: The [Seurat manual](https://cran.r-project.org/web/packages/Seurat/Seurat.pdf) does a good job explaining the parameters for any function.

> **min.cells** Include features detected in at least this many cells. Will subset the counts matrix as well. To reintroduce excluded
> features, create a new object with a lower cutoff.
>
> **min.features** Include cells where at least this many features are detected.

The values they picked here are somewhat arbitrary, but `min.cells` helps limit the number of genes used by removing those unlikely to play any part in differentiating groups of cells due to being expressed in very few cells. In general, most genes removed will be those with zero counts across all cells. `min.features` removes dead cells cells and empty droplets where few genes are detected.
Question: Hi all,

I have BioNano Scaffolds of sequence data (no solid reference at the moment) and I am trying to orient them in the correct way to so I can do PCR and attempt to complete the sequence. One of my scaffolds is 67 MB, but I believe (based on BLATs and other context clues from the project) that this sequence is actually reversed. Is there a way using samtools or the command line to easily flip this sequence? Given that it is so big I would rather stay on the cluster and not have to download it to my computer to try and flip - but if that is what I need to do so be it.

Thanks so much!

Answer: If you have the emboss package installed on the server you can use `revseq seq.fa -reverse -nocomplement -outseq seq.rev.fa`.
Question: hi everyone,

I have just started learning genomics as a part of my bioinformatics degree and I've been introduced to linux for handling fastq files and using fast qc. can anybody suggest some good learning resources which is more inclined towards linux for genomics.

Answer: This [**UNIX primer**][1] from Korf lab is an excellent resource to get started.

Titus Brown's [**NGS data analysis workshop**][2] is also a great resource.

[1]: http://korflab.ucdavis.edu/Unix_and_Perl/current.html#part1
[2]: https://angus.readthedocs.io/en/2019/toc.html
Question: Hi all.

Recently, I got sequence for H3K27ac ChIP-seq data performed in various stimulation and knockdown condition. I have processed the reads and did peak calling using macs2. Looking into IGV browser, the enrichment looks good and the peaks are identified. The peaks called are of in the size range from ~1kb to 10kb. My concern here is whether the size of peaks identifies is really broad or is it ok to get large size peaks from H3k27ac chip. Next if want to perform motif identification, should I take all the broad peaks or perform peak calling to get broad peaks of lesser size (~1kb).

I would highly appreciate your suggestions and feedbacks

Thanks

peak calling method:

ls bowtie_out/*.25M.bam | cut -d'/' -f2 | parallel --verbose 'macs2 callpeak -t bowtie_out/{} -c ../Input_R1.bam -g mm -n {=s/.25M.bam//g=} --verbose 2 --broad --broad-cutoff 0.01 --cutoff-analysis --fe-cutoff 4 --outdir macs_peak'

IGV browser snapshot:

Zoom of portion highlighted in above image:
free image hosting

Answer: If the signals in IGV are looking reasonable over background for ~10kb domains, it is totally acceptable. There are something called **Super-enhancers** which are minimum of 10-12kb in size with high H3K27ac signal. Just google the term, you'll find number of papers solely based on super-enhancers.

Regarding motif analysis, instead of using total peak length, identify nucleosome free region within your H3K27ac peaks. These are potentially identified by small valleys which lack signal within peaks. `Homer` tool has a function called `nfr` which is what you want.

![](http://homer.ucsd.edu/homer/ngs/peaks.nfrExample.png)

BTW, is there a specific reason to use `broad` option with MACS in your case? In my experience, narrowpeaks from MACS are of good quality and works very well for all types of analyses.
Question: This is not a very memory friendly way of counting sequences from a multi fasta, any ideas to improve this?

generator = SeqIO.parse("test_fasta.fasta","fasta")
sizes = [len(rec) for rec in SeqIO.parse("test_fasta.fasta", "fasta")]

I'm avoiding using tools like grep since I want to make this more portable

Answer: Standard Python will be faster than BioPython:

fh = open("test_fasta.fasta")
n = 0
for line in fh:
if line.startswith(">"):
n += 1
fh.close()

or shorter and possibly faster:

num = len([1 for line in open("test_fasta.fasta") if line.startswith(">")])
Question: Hi Biostars

I am trying to convert sra from [PRJNA282735][1] dataset to fastq and I am getting following error...

> fastq-dump.2.1.7 fatal: SIGNAL - Segmentation fault

My fastq-dump command is

fastq-dump --split-3 SRR2016445.sra -O SRR2016445

I am not able to find similar error elsewhere. The [ENA page][2] for some samples of this dataset has three files per SRX experiment (e.g. SRR2016445.fastq, SRR2016445_1.fastq and SRR2016445_2.fastq).

This is unusual for me as I usually get one or two SRR runs per experiment (depending on single end paired end) but never 3. I am wondering if this is the reason for getting errors.

Anybody with similar experience?

[1]: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA282735/
[2]: http://www.ebi.ac.uk/ena/data/view/PRJNA282735

Answer: I think 2.1. is fairly old, with a slightly newer version 2.4.2 I get:
fastq-dump.2.4.2 err: error unexpected while resolving tree within virtual file system module - failed to resolve accession 'SRR2016445' - Obsolete software. See https://github.com/ncbi/sra-tools/wiki ( 406 )

The latest release on github is 2.8, you should get and install the latest version as outlined here: https://github.com/ncbi/sra-tools/wiki
Question: Dear All,

I would like to ask you about an issue am facing. While doing exome sequencing I earlier did an estimation to what extent my reads were falling on the target regions, ie the target intervals that are used for target enrichment and I found for my samples the reads were around 75% spanning the exonic region. But now when am translating them into the variants and trying to find the somatic variants and annotating the somatic variants with all 3 annovar, oncotator and snpEff I find only 30% SNPs(novel ones which are not in dbSNP) are actually on the exons. The rest are spanning the intronic,intergenic,splice ,UTRs etc. Is this a likely scenario? How often do you see the SNPs mostly annotated in the non-exonic regions even when you use the exome sequencing provided your reads have high coverage in the exonic intervals used for target enrichment? I would like to know your advice in such cases.

Thanks

Answer: Does this include synonymous mutations? If you filtered out synonymous exonic mutations you'll find that is a large contributor as well. Most of the non-exonic regions you listed are all either targeted, or will be captured due to overlap (some areas extending into introns) and the majority of all mutations in a cell happen to occur in non exonic regions simply because exonic regions make up such a small percentage of the genome. Of course since you have filtered out off-target variants in this case, you will enrich for exonic variants.

Coupled with filtering out mutations that are also in dbSNP (which many valid somatic mutations may overlap or match identically mutations found in dbSNP by the way, plenty of known somatic variants that appear in COSMIC, TCGA, etc are also found in dbSNP) the number you are reporting (30%) doesn't surprise me that much.
Question: I am trying to calculate the coverage of a ChIP data set over a very limited set of gene coordinates. I have very simple bed files for both (after trying to eliminate any cause of the error) that are just the bare bones "chr start end" columns, and every time I run the command:

bedtools coverage -a a.bed -b b.bed

I get this error:

```
ERROR: Received illegal bin number 4294967295 from getBin call.
ERROR: Unable to add record to tree
```

I used grep to search for this number and it's nowhere to be found in either of my files, which makes sense given that my genome isn't that big. So why am I still getting this error?

Thanks for your input!

Answer: The BED files look reasonable.

Maybe there are invisible characters somewhere in the file? Can you try running `mac2unix` or `dos2unix` on them?
Question: Hi all,

I was just wondering if anybody had any experience with coloring something like a UMAP made in ggplot based on the expression of multiple genes at the same time? What I want to do is something like the blend function in Seurat featureplots, but with 3 genes / colors instead of 2.

I'm looking to make something like this:

![Image 1][1]

Where the colors for the genes combine where there is overlap.

What I've gotten to so far is

ggplot(FD, vars = c("UMAP_1", "UMAP_2", "FOSL2", "JUNB", "HES1"), aes(x = UMAP_1, y = UMAP_2, colour = FOSL2)) +
geom_point(size=0.3, alpha=1) +
scale_colour_gradientn(colours = c("lightgrey", colour1), limits = c(0, 0.3), oob = scales::squish) +
new_scale_color() +
geom_point(aes(colour = JUNB), size=0.3, alpha=0.7) +
scale_colour_gradientn(colours = c("lightgrey", colour2), limits = c(0.1, 0.2), oob = scales::squish) +
new_scale_color() +
geom_point(aes(colour = HES1), size=0.3, alpha=0.1) +
scale_colour_gradientn(colours = c("lightgrey", colour3), limits = c(0, 0.3), oob = scales::squish)

Where FD is a data frame containing the information from the seurat object for the UMAP coordinates and the expression levels of the three genes of interest. All I can get is a plot where the points from one layer obscure those below it, I've tried messing around with the colours, gradients, alpha and scales but I'm guessing I'm doing it the wrong way.

![Image 2][2]

If anyone knows of a way to make this work or has any suggestions on something else to try that would be very much appreciated.

[1]: /media/images/e6f45dd8-da0d-4858-950c-74d3f4c4
[2]: /media/images/a661c7cc-7b77-4efe-b921-e1eb4412

Answer: Just in case anyone was searching for an answer for this in the future, I asked the same question over on StackExchange too and was suggested an answer that worked perfectly for me.

https://stackoverflow.com/questions/69609056/colour-umap-based-on-expression-of-multiple-genes-in-ggplot2/