cannin/biostars_qa huggingface dataset:
Question: Hi guys. I'm working with a data set and I've written a for loop that generates barplots for my input tRNA data. I've been able to generate my bar plots. However, I have a column with my P values and I want to annotate on my bar plots any significant log fold changes (Padj<=0.05) with an asterisk, and insignificant values (Padj>0.05) with NS. How do I do this and incorporate it into my for loop. Been stuck on this all weekend. THANKS IN ADVANCE!!!

***head(df) ...

Isodecoder Anticodon Loci Fragment_type Label LFC Padj
1 Ala AGC 1 wholecounts Ala-AGC-1-wholecounts 3.6111408 8.890000e-10
2 Ala AGC 10 wholecounts Ala-AGC-10-wholecounts 0.3910508 6.404589e-01
3 Ala AGC 12 wholecounts Ala-AGC-12-wholecounts NA NA
4 Ala AGC 2 wholecounts Ala-AGC-2-wholecounts 0.2785932 4.741744e-01
5 Ala AGC 3 wholecounts Ala-AGC-3-wholecounts 1.3492899 3.269910e-04
6 Ala AGC 4 wholecounts Ala-AGC-4-wholecounts -0.4200119 1.475002e-01

***This is my current code:

df <- read.delim("/mnt/data/wholecounts_barplot_input.txt", stringsAsFactors = FALSE)

library(ggplot2)

plot_list <- list()

for (iso in c("Ala", "Arg", "Asn", "Asp", "Cys", "Gln", "Glu", "Gly", "His", "Ile", "Leu", "Lys", "Met", "Phe", "Pro", "SeC", "Ser", "Thr", "Trp", "Tyr", "Val", "iMet")) {

df2 <- subset(df, Isodecoder == iso)

temp_plot <- ggplot(data = df2, aes(x = Label, y = LFC, fill = Anticodon))+
geom_bar(stat = "identity")+
xlab("tRNA")+
ylab("Log2FC")+
ylim(c(-2, 4))+
theme_bw()+
theme(axis.text.x = element_text(size = 12, angle = 45, hjust = 1), plot.margin = margin(0.5,0.5,0.5,2, "cm"))+
ggtitle(iso, "Full Length tRNA LFCs")

plot_list[[iso]] <- temp_plot

pdf(paste(iso, "wholecounts_LFCs.pdf", sep = "_"), bg = "white", width = 15, height = 7)
print(temp_plot)
dev.off()
}

Answer: Annotate bars with *geom_text*:

ggplot(data = df1, aes(x = Label, y = LFC,
fill = Anticodon,
label = ifelse(Padj < 0.05, "*", "NS")) # <--- See here
) +
geom_bar(stat = "identity") +
geom_text(vjust = 0) + # <----------------------------------------- and here
xlab("tRNA") +
ylab("Log2FC") +
ylim(c(-2, 4)) +
theme_bw() +
theme(axis.text.x = element_text(size = 12, angle = 45, hjust = 1),
plot.margin = margin(0.5, 0.5, 0.5, 2, "cm"))
Question: Hi all,

any advice on where i can find the genetic map for the GRCh38 version of 1000Genomes? this is needed as a parameter in the imputation process by IMPUTE2.

Answer: hi all,

i finally found the genetic map for build GRCh38. it is at the below link if anyone needs to have it

https://alkesgroup.broadinstitute.org/Eagle/downloads/tables/
Question: Hi,

I have a read counts file from single cell seq from RSEM, I want to store that in a matrix.mtx format required for Seurat pipeline, how can I do that?

> head(all.counts)
HSC1_1 HSC1_2 HSC1_3 HSC1_4 HSC1_5 HSC1_6 HSC1_7
0610005C13Rik 0 0 0 0 0 0 0
0610007P14Rik 25 0 304 0 154 550 117
0610008F07Rik 0 0 0 0 0 0 0
0610009B14Rik 0 0 0 0 0 0 0
0610009B22Rik 0 0 0 8 83 521 236
0610009D07Rik 241 0 66 60 356 173 385

to

%%MatrixMarket matrix coordinate real general
%
32738 2700 2286884
32709 1 4
32707 1 1
32706 1 10
32704 1 1
32703 1 5
32702 1 6
32700 1 10
32699 1 25
32698 1 3
32697 1 8
32527 1 1

Answer: Hi!
Use the Matrix library to create and write a sparse matrix:

library(Matrix)

# generate single-cell RNA seq data
gbm <- replicate(10, rpois(10, 1))
rownames(gbm) <- paste0("gene_", 1:nrow(gbm))
colnames(gbm) <- paste0("cell_",1:ncol(gbm))

# save sparse matrix
sparse.gbm <- Matrix(gbm , sparse = T )
head(sparse.gbm)
writeMM(obj = sparse.gbm, file="matrix.mtx")

# save genes and cells names
write(x = rownames(gbm), file = "genes.tsv")
write(x = colnames(gbm), file = "barcodes.tsv")

Result:

$ head matrix.mtx
%%MatrixMarket matrix coordinate integer general
10 10 65
1 1 3
2 1 1
3 1 3
6 1 1
7 1 1
8 1 1
10 1 3
1 2 2

Jean-Baptiste
Question: Hi,

I have a set of samples and log10 total read counts of them in a data frame like below

> head(a)
log10_total_counts
A1 6.468503
A10 6.565213
A11 6.752139
A12 5.078598
A2 6.277342
A3 6.473411
>

![enter image description here][1]

[1]: https://i.ibb.co/P6mmKc5/Picture1.png

For instance red dots have less than 1000000 counts and black dots have more thsn 1000000 counts.

How could I repreduce this plot please?

Answer: Use `ggplot2`.

Something is missing in your example; the labels at the bottom do not match you data frame, and it's unclear how your sample dataframe would count up to values near 1000000. You might want to precompute/add a QC vector based on that to your dataframe; basically have your dataframe in a state that's ready to be plotted with minimal calculations.

Here's an example; I generate a dataframe with normally distributed values centered on 6.0 (seemed close to what your example had). I plot a dotted line at 6; and anything below 6 becomes red.

![Example ggplot call with geom_points and two colours. ][1]

For the style, I directly took the function `theme_Publication()` from this [link][2].

Here's the code:

library("ggplot2")

# Themes for style
theme_Publication <- function(base_size=14, base_family="helvetica") { ... }
scale_fill_Publication <- function(...){ ... }
scale_colour_Publication <- function(...){ ... }

# Function to plot points as described
plotReads <-
function(DF) {
ggplot(DF) +
theme_Publication() + # Style it
geom_point(mapping = aes(x = as.factor(1:nrow(DF)), # Sequentially
y = log10_total_counts, # Y is the count values
color = QC)) + # Color it by the QC vector in the dataframe
scale_color_manual(values = c("PASS" = "black", "FAIL" = "red")) + # Tell it what colours to use.
geom_hline(yintercept = 6.0, linetype="dotted") + # Add the dotted line at 6.0
xlab("Sample") + # Label the X
ylab("log10(Total Reads)") + # Label the Y
ylim(c(-0.2, 8.2)) # Extend the Y axis to be between this range.
}

DF = data.frame(
log10_total_counts = rnorm(30, mean = 6, sd = 0.5), # Randomized data
QC = as.factor(c( rep("FAIL",15), rep("PASS", 15) )) # A vector that is just half FAIL and half PASS
)
plotReads(DF)

DF$QC <- as.factor(ifelse(DF$log10_total_counts < 6.0, "FAIL", "PASS")) # Now actually add some logic to what should be PASS or FAIL.
plotReads(DF)

[1]: https://i.imgur.com/zqbjj7a.png
[2]: https://rpubs.com/Koundy/71792
Question: Hi,

I ran CRISP to do the variant calling in some Drosophila pooled samples and when I try to filter for minimum allele frequency with vcftools I get the following error:

VCFtools - 0.1.14
(C) Adam Auton and Anthony Marcketta 2009

Parameters as interpreted:
--vcf new_variantcalls.vcf
--maf 0.15
--out berg

After filtering, kept 14 out of 14 Individuals
Error: Polyploidy found, and not supported by vcftools: YHet:49

Does anybody know how to fix this?

I run the grep command Pierre suggested and this is the output is here: https://jpst.it/MUgG

Answer: looking at the vcftools source code: https://github.com/kiwiroy/vcftools/blob/master/cpp/vcf_entry_setters.cpp#L97

if (pos != pos2)
LOG.error("Polyploidy found, and not supported by vcftools: " + CHROM + ":" + header::int2str(POS));

you have a strange genotype in CHROM= YHet at POS = 49. (something that doesn't look like '([0-9]+/[0-9]+)'.
I suppose you have more than two alleles for a Genotype

please, show us the output of

grep -v "^#" new_variantcalls.vcf | cut -f 10- | cut -d ':' -f 1 | sort | uniq

> Does anybody know how to fix this?

change the way your VCF is called. ...
Question: My aim is to create a bash loop over my fastq files that are in a directory like so:
https://ibb.co/cMvx7a
whereby, paired read files are organised in the form:

1

1

2

2

3

3

etc..

I don't know a way of dealing with these paired files within the cutadapt loop
Something along these lines is what I'm looking for:

for file in folder

do
cutadapt -a ADAPTER_FWD -A ADAPTER_REV -o out.1.fastq -p out.2.fastq reads.1.fastq reads.2.fastq
done

but I want to do pairs together from the folder and I don't know how.

Thanks in advance.

Answer: This is not a bioinformatics question, it is more a shell scripting question.

This draft should help you to get started:

for i in *_R1.fastq.gz
do
SAMPLE=$(echo ${i} | sed "s/_R1\.fastq\.gz//")
echo ${SAMPLE}_R1.fastq.gz ${SAMPLE}_R2.fastq.gz
done
Question: While trying to run RepeatMasker I get the error

The RepeatMasker.lib file is out of date. This version of RepeatMasker requires library version 20140131 or higher. Your version is 20110419.

I am still waiting for a response from [Genetic Information Research Institute][1] for a password to login for downloading and using their libraries to run repeatmasker. Would that fix this error?

Also I haven't been able to see the default libraries but I am guessing that it is suitable to be used on any genome when we haven't got our own custom made libraries. Is that so?

[1]: http://www.girinst.org/

Answer: After I received my password from giri I downloaded the new repeatmaster file, unzipped it into the libraries folder that repeatmaster had built and I have had no problems using it since then. I have not however tried to install a custom library yet.

Susanne
Question: I have a bam file produces by BWA-MEM. It has some unmapped reads, e.g.

D00733:389:CD1T7ANXX:3:1101:1572:2235	77	*	0	0	*	*	0	0	CAGTTTCACTGTATAAATTGCTTATACTTAGACATGCATGGCTTAATCTT	AAB=AFGDGCFFGGGGGGGGGCGGGGGGGGGGGGGGGGGGGGGGGGGGGF	AS:i:0	XS:i:0
D00733:389:CD1T7ANXX:3:1101:1572:2235	141	*	0	0	*	*	0	0	GTATCTTCTAGAGAGAGGGAATGGGCGAGAGAAAAAGAGATTTCGGTTTC	BBB@BGGGGGGGGFGGGGGGGGEGGGGGGDGFGGGGGGGGEGGGGGFGGG	AS:i:0	XS:i:0
D00733:389:CD1T7ANXX:3:1101:6797:2243	77	*	0	0	*	*	0	0	TGTCTGGACCTGGTGAGTTTCCCCGTGTTGAGTCAAATTAAGCCGCAGGC	3A<0BDGGGGGGGGGGGGGFGGGGGFGGGGGGGGGGGGGGGGGGGGGGGG	AS:i:0	XS:i:0

I use pysam to count some stats on the bam file, but for some reason pysam does not find these unmapped reads.

bam=pysam.AlignmentFile(file,"rb")
for line in bam.fetch():
line=line.tostring(bam)
line=line.split("\t")
if line[2]=="*":
print(line)

The code does not return anything

Any ideas how to fix this?

Thanks

Answer: you should add `until_eof=True` in the `fetch()`
Question: I am using Picard [MarkDuplicates][1] to mark duplicates in my position sorted BAM. Then, I variant call using bcftools [mpileup][2] and bcftools [call][3] like so:

bcftools mpileup -f reference.fa positionsort.bam | bcftools call -mv -Ob -o calls.bcf

To my knowledge, bcftools [mpileup][2] ignores duplicates by default ([here][4] and [here][5]). Initially I assumed that this would be the same case for bcftools [call][3]. However, I was skeptical and decided to do a test where I variant call when I mark duplicates with Picard, and where I variant call without marking duplicates. To my surprise, this resulted in both cases in the exact same number of variants. So, does this mean that bcftools call *doesn't* ignore duplicate markings? Or have I made a mistake?

**edit:** I did in fact make a mistake, see comments below

[1]: https://gatk.broadinstitute.org/hc/en-us/articles/360037052812-MarkDuplicates-Picard-
[2]: https://samtools.github.io/bcftools/bcftools.html#mpileup
[3]: https://samtools.github.io/bcftools/bcftools.html#call
[4]: https://www.biostars.org/p/3917/
[5]: https://www.biostars.org/p/18182/

Answer: So I did find a bug in my pipeline that caused the variant calling to be done only on the marked duplicates files, and not on the files where I didn't mark the duplicates, explaining why I found the exact same number of variants in both of these cases. I fixed that mistake, did some more tests and also tested SAMtools [markdup][1]. This resulted in the following:
![table][2]

I tested marking the duplicates, and removing them, and then doing the variant calling. The number of SNPs called with BCFtools are shown in the table above. As you can see, the # of variants called when marking duplicates and removing duplicates are basically the exact same, except for Picard SAMPLE5, where removedup resulted in 1 variant less then markdup (in red in the table). So, I think it's safe to say that bcftools [call][3] also ignores duplicate markings, just like [mpileup][4].

[1]: http://www.htslib.org/doc/samtools-markdup.html
[2]: /media/images/890f6bc8-b449-41cb-b44f-5b6d5919
[3]: https://samtools.github.io/bcftools/bcftools.html#call
[4]: https://samtools.github.io/bcftools/bcftools.html#mpileup
Question: I am using [snp.plotter][1] to illustrate LD heatmaps along side a zoom plot. [Here][2] is an example of what I mean.

However, in [this manuscript][3], they have images that look to be made by snp.plotter, but have gene diagrams overlaid, like can be done in LocusZoom. In [the documentation for snp.plotter][4], however, I cannot find any option to add such a diagram.

Does anyone know how this is being accomplished? I know how to use Image Magick if an ad hoc addition is required, but at present I am not sure if the diagram is being added through the program or by another means...

[1]: https://github.com/cannin/snp_plotter
[2]: https://camo.githubusercontent.com/6c7d77ead2dbe3a360124ac52ff12422c1156b3e/68747470733a2f2f7261772e6769746875622e636f6d2f63616e6e696e2f736e705f706c6f747465722f6d61737465722f646f63732f6e756c6c5f7273715f706879735f636f6c6f722e706e67
[3]: http://www.nature.com/ng/journal/v46/n9/fig_tab/ng.3052_F1.html
[4]: https://cran.r-project.org/web/packages/snp.plotter/snp.plotter.pdf

Answer: You're correct. snp.plotter does not include a feature to incorporate the LocusZoom plots. Likely, the authors of the article created the two plots and then combined them in Illustrator or Inkscape. snp.plotter can produce PDFs so you would not lose quality by combining the images.
Question: I'd like to generate a plot similar to the included image. What is the name of such a plot and does there exist a simple tool to generate, preferably compatible with multiqc and/or R/bioconductor?

![enter image description here][1]

[1]: https://external-content.duckduckgo.com/iu/?u=https%3A%2F%2Fwww.takarabio.com%2FImages%2F150-cDNA_Synthesis%2FCC-TechFigs%2FSINGLE_TN04_%2520109256.jpeg&f=1&nofb=1

Answer: Maybe [CoverageView][1]? Looks like it was made for other experiment types but since it uses bam files as input it might work. Else Bioconductors bamsignals have a [bamCoverage()][2] function

Else I know RSeQC have a gene-coverage tool which you can find [here][3].

Just be aware that generally in single cell data (especially for 5' end data such as 10X) you have way to few reads per gene to look at coverage. Most genes have very few reads so how much sense it makes to calculate this I'm note sure.

[1]: https://bioconductor.org/packages/release/bioc/html/CoverageView.html
[2]: https://www.bioconductor.org/packages/release/bioc/vignettes/bamsignals/inst/doc/bamsignals.html#read-coverage-with-bamcoverage
[3]: http://rseqc.sourceforge.net/#genebody-coverage-py
Question: Hi all,

I was recently given a data set (run through Illumina in 2014) to merge with a set run through Illumina this year. The set I just received was not demultiplexed and was dual barcoded. I know that Illumina's bcl2fastq can handle dual barcoded sets and demultiplex, but I don't have any of the BaseCall data it typically uses to demultiplex. I only have .fastq.gz files and a mapping file to work with. Does anyone know if I can still use the bcl2fastq without the BaseCall data, or if there is something that would work better?

tldr; I need to demultiplex dual barcoded .fastq.gz files

Answer: You can't use bcl2fastq for this.

The bad news is that there might not be a tool to do as it is a task that is usually handled at the instrument level so there is less of a need to do it.

You would probably need to use a tool designed for cutting adapters to filter the reads:

http://cutadapt.readthedocs.io/en/stable/guide.html#filtering-reads

if it is a single end sequencing then this might be sufficient. If the reads are paired-end you'd probably need to process/reconcile both pairs to keep them in sync.
Question: Hello,

For more than a month I am trying many tools for finding DE miRNAs in fungi but error, error, error

online tools don't support fungi and local tools were mostly written in perl needed mirdeep2 algorithms

believe my failure is not because of my laziness, I tried whatever I could

although being too rude with a great expectation, please someone help me to install **mirdeep** via TeamViewer or whatever

I asked institute admin, he tried his best but my problem was not solved

thank you in advance

Answer: Hi Angel,
I dont know if there is any specific reason that you are using miRDeep2. But miRDeep* claims that it is better (in predicting novel miRNAs) than miRDeep2 (http://www.australianprostatecentre.org/research/software/mirdeep-star) and has GUI. You can try this.
Additionally, there is a web-server oasis2.0 (https://oasis.dzne.de/index.php) that can be used for miRNA identification, DE study and Target prediction. Though your organism of interest is not listed there, there is an option "All Genomes of miRBase v21" under Reference Genome drop-down. I am not sure, if you can use this option.
Question: I am analyzing rna seq data from illumina stranded protocol 126bp -PE , sequences have nextera adapters using cutadapt I trimmed all the sequences but now I have reads with different lengths, also my insert size form picard has weird peaks, find below has any one experienced this before?
![enter image description here][1]

#EDIT

Apologies I did little more digging around the pipeline and found out that flexbar was used to remove adapters

and the command used was

flexbar --adapters adapters.file --adapter-trim-end RIGHT --length-dist --threads 12 --adapter-min-overlap 7 --max-uncalled 250 --min-read-length 25

adapter used
>Read_1_Sequencing_Primer_3_to_5
AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT
>Read_2_Sequencing_Primer_3_to_5
AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC

#Read1
![enter image description here][2]

#Read2![enter image description here][3]

#Using leeHom tool.

![enter image description here][4]

#log file leeHom

Total reads :231670526 100%
Merged (trimming) 22238730 9.59929%
Merged (overlap) 0 0%
Kept PE/SR 93263449 40.2569%
Trimmed SR 0 0%
Adapter dimers/chimeras 533995 0.230498%
Failed Key 0 0%

[1]: https://i.imgur.com/224uGAh.png
[2]: https://i.imgur.com/cMBpvVe.png
[3]: https://i.imgur.com/2SmFqfB.png
[4]: https://i.imgur.com/y4Vl20a.png

Answer: If you try leeHom ([https://grenaud.github.io/leeHom/][1]), it can trim adapters and merge overlapping using a Bayesian algorithm so no cutoffs are required. It will produce 3 types of files to map:

1. [PREFIX].fq.gz	Sequences that were trimmed and merged confidently by leeHom
2. [PREFIX]_r1.fq.gz	Forward reads that were neither trimmed nor merged confidently by leeHom
3. [PREFIX]_r2.fq.gz	Reverse reads that were neither trimmed nor merged confidently by leeHo

Map the .fq.gz as single reads and r1/r2 as pairs. To plot the insert size, get pysam and run the following script as such: python insertSize.py in.bam

#!/usr/bin/python

import pysam
import sys;

bamInputFile = pysam.Samfile(sys.argv[1], "rb");

for read in bamInputFile:
if read.is_unmapped :
continue;
if read.is_paired :
if read.is_proper_pair and read.is_read1 :
print abs(read.template_length);
else:
print len(read.query_alignment_sequence);

bamInputFile.close();

[1]: https://grenaud.github.io/leeHom/
Question: Is there a difference between the nuccore and nucleotide entrez databases? In my simple tests they seem to return the same data. However [the biopython tutorial][1] makes use of both of them.

My testing from python:

```py
Entrez.read(Entrez.einfo(db='nuccore'))["DbInfo"]["Description"]
'Core Nucleotide db'

Entrez.read(Entrez.einfo(db='nuccore'))["DbInfo"]["Count"]
'202831776'

Entrez.read(Entrez.einfo(db='nucleotide'))["DbInfo"]["Description"]
'Core Nucleotide db'

Entrez.read(Entrez.einfo(db='nucleotide'))["DbInfo"]["Count"]
'202831776'
```

I also tried using efetch and didn't see any differences in the returned data.

[1]: http://biopython.org/DIST/docs/tutorial/Tutorial.html

Answer: When trying to figure out the difference between the two databases, I also ran the same commands on the database descriptions, which indicates to me that they are the same database.

In the Biopython tutorial under Chapter 9, it mentions "The `Bio.Entrez` module makes use of the Entrez Programming Utilities (also known as EUtils)." According to NCBI's book "Entrez Programming Utilities Help", the nuccore database is just E-utility database name corresponding to [the Nucleotide Entrez database](http://www.ncbi.nlm.nih.gov/books/NBK25497/table/chapter2.T._entrez_unique_identifiers_ui/?report=objectonly).

Additionally, going to the Nucleotide database on NCBI, I notice "nuccore" in [the web address](http://www.ncbi.nlm.nih.gov/nuccore).

So I think they're the same database :)
Question: This seems simple maths behind, But I always get confused with the terminology of **"Fold change", "Log Fold Change" *(LFC)*, and "Percentage of Change"**., When I see some of the figure descriptions given in **volcano plots**, specially X axis, I get confused.

i.e. [Page 64][1]

When I google, I found [this figure][2] is easy to understand, but the figure description does not make sense for me.

for e.g.
If we take two-fold change;
isn't this the way we calculate LFC;

LFC = Log2 (Foldchange) = Log2(2) =>1

if X axix is in Log2 fold change, shouldn't two-fold change denotes at **± 1**? Why they denote 2 fold change cutoff line in **± 2**?

-----------------
**Below VP; Fold Change greater than 2 times and significant P value less than 0.05**

See this example:

![enter image description here][3]

and this:

![enter image description here][4]

Could someone help me to understand this terminology ( **"Fold change", "Log Fold Change" *(LFC)*, and "Percentage of Change"**. ) easily with an example?

[1]: http://www.nathalievialaneix.eu/doc/pdf/tutorial-rnaseq.pdf
[2]: https://documentation.partek.com/display/FLOWDOC/Volcano+Plot
[3]: https://cdn.xlstat.com/uploads/assets/Tutorials/DifferentialExpression/differentialExpression06.png
[4]: http://chenyuan.date/images/Volcano-plot-unnamed-chunk-4-1.png

Answer: Hey @ Wuschel

These answers here helped me a lot in the beginning to understand these terminologies that I even saved in my browser bookmarks.

[Question: Difference between FC and LogFC][1]

[Dumb question about LogFC][2]

[Question: Report and compare log2FC vs FC][3]

I think it can help you about log2FC and FC, but about "% of change" I never read anything about it.

Another thing, if you want to learn how to create beautiful volcano plots read this tutorial "[Tool: EnhancedVolcano: Publication-ready volcano plots with enhanced colouring and labeling][4]"

[1]: https://www.biostars.org/p/312980/
[2]: https://www.reddit.com/r/labrats/comments/7odtki/dumb_question_about_logfc/
[3]: https://support.bioconductor.org/p/66163/
[4]: https://www.biostars.org/p/335751/
Question: I am trying to understand Gene expression data analysis. Could someone the percentage of change of the GE calculating from log2 fold change (L2FC)?

> FoldChange = μ_Treatment / μ_Control
>
> log2FC=Log2(μ_Treatment)-Log2(μ_Control)

Answer: [calculate fold change (FC)][1]
When comparing these log transformed values, we use the quotient rule of logarithms:
>log(A/B) = log(A) - log(B)

>log(A) = 4

>log(B) = 1

>Therefore:

>log(A/B) = 4 - 1

>log(A/B) = 3 This gives a 3-fold change.

Please note that in this case we are reporting the log(fold change). Biologists often use the log(fold change) because without taking the log, down regulated genes would have values between 0 and 1, whereas up regulated genes would have any value between 1 and infinity. This distribution makes graphing and further statistical analysis difficult. Taking the log typically makes the resulting values more normally distributed, which is better for further analysis.

[1]: https://ucsc-xena.gitbook.io/project/faq/data-and-datasets#how-do-i-calculate-fold-change-fc