cannin/biostars_qa huggingface dataset:
Question: Example input: multi-sample VCF (adapted from [www.internationalgenome.org][1]):

**Note:** my actual file is bgzipped and tabixed with ~2Mln variants (rows) and ~1000 samples (columns).

##fileformat=VCFv4.0
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=1000GenomesPilot-NCBI36
##phasing=partial
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##FILTER=
##FILTER=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA00003
20 14370 rs6054257 G A 29 PASS NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ 0/0:48:1:51,51 1/0:48:8:51,51 1/1:43:5:.,.
20 17330 . T A 3 q10 NS=3;DP=11;AF=0.017 GT:GQ:DP:HQ 0/0:49:3:58,50 0/1:3:5:65,3 0/0:41:3
20 1110696 rs6040355 A G,T 67 PASS NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQ:DP:HQ 1/2:21:6:23,27 2/1:2:0:18,2 2/2:35:4
21 1230237 . T . 47 PASS NS=3;DP=13;AA=T GT:GQ:DP:HQ 0/0:54:7:56,60 0/0:48:4:51,51 0/0:61:2
21 1234567 microsat1 GTCT G,GTACT 50 PASS NS=3;DP=9;AA=G GT:GQ:DP 0/1:35:4 0/2:17:2 1/1:40:3

Expected output:

**20.txt**

NA00001	NA00002	NA00003
0/0	1/0	1/1
0/0	0/1	0/0
1/2	1/2	2/2

**21.txt**

NA00001	NA00002	NA00003
0/0	0/0	0/0
0/1	0/2	1/1

I was thinking of using `cut | sed` combo with some `regex`, but thought there must be already some tool out there, maybe *[bcftools][2]* (couldn't get the right flags to work) ?

Any other ideas?

[1]: http://www.internationalgenome.org/wiki/Analysis/Variant%20Call%20Format/vcf-variant-call-format-version-40/
[2]: http://www.htslib.org/doc/bcftools.html

Answer: From your expected output, it is evident that you wish to not split multi-allelic records. I strongly recommend splitting them as tools become more predictable with normalized variants.

You can use

bcftools query -r -Hf '[ %GT]\n' # to get genotypes

to get this output, but beware, the header will be a bit wonky, in the format `[]:GT`. You can `sed` replace the `\[\d+\]([^:]+):GT` with `\1` to clean up the header and get to your desired output.

for chr in 20 21
do
bcftools query -Ov -r ${chr} -Hf '[ %GT]\n' | sed -r '1s#\[\d+\]([^:]+):GT#\1#g' >${chr}.txt #command untested, use with caution
done
Question: Hi, I used GATK germline variant calling pipeline to call short variants on paired end fastq files. After got the final analysis ready vcf, applied some extra filters, I inspected bam files in IGV for those variants of interest and found some strange things for one sample. Two variants of interest in this sample can only be found on inversion reads.

In this first graph, the alternative allele G can only be found in RR and LL reads (blue color) in IGV. 13 out of 15 inversion reads have this G allele.
![enter image description here][1]

In the second graph, similarly, the alternative allele T can only be found in inversion reads. All the inversion reads have this T allele.
![enter image description here][2]

Further, I realized that all the inversion reads have same size shown in figure 3.
![enter image description here][3]

I wonder if these inversions are true inversions or they are artifacts (given all the same size) thus the variants only found on these reads are also not real.

[1]: /media/images/33ddd7ec-6365-4126-859f-b9b9b990
[2]: /media/images/fbff412d-756e-4f82-b0fc-72be07ff
[3]: /media/images/623bd4b4-df4e-4fad-a59e-fb8f8c92

Answer: Hi xukeren,

Not every LL or RR read shows a true inversion. I think you have to further inspect your reads in IGV if you want to understand what is going on.

It could be a ambiguous mapping due to the short length of the reads. What is the mapping quality of the reads? You should see it if you click on a read.

You should activate seeing mismatched bases if you have soft clipped bases: View->Preferences->Alignments->Show mismatched bases. Your short reads might get longer tails of colorfull sequences which were not used in mapping at this position, but might have been used for a supplementary alignment at another position. (You might have to reload your BAM to IGV to take effect)

You need to find out where the corresponding reads for each of your RR/LL reads are:
Click on a specific read and look at the information in the section about "Mate".
To easily visualize it, you can try turning on "View as pairs" in the rightclick menu of the BAM-file. If the corresponding reads are close enough they get connected by slim lines.
The other way would be "Go to mate" after right-clicking on a specific read (only available if "View as pairs is turned off"). IGV will jump to the mapping position of the corresponding read and highlight both in a unique color.
Question: Hi,

I'm trying to pipe the following commands for ChipSeq analysis:

${bowtie2_source} -x ${ref_genome} -U ${fastq_file} -S - | ${samtools} view -Shub - | ${samtools} sort - | ${samtools} rmdup -S - - ${target_dir}/${sample_name}.bam

Up till the sorting part it works fine, but in the end, no bam file is created and it gives me the message:

[bam_rmdup] input SAM does not have header. Abort!

What is the problem? and How should I fix it?

Thanks!

Answer: In the last step you seem to specify two `-` flags, I think that might attempt to open the stream twice, and since the second time it is empty, there is no header information.

${samtools} rmdup -S - - ${target_dir}/${sample_name}.bam

In general, in such cases decompose the command into individual steps, each operating on a differently named file, then investigate the resulting file after each step.

The samtools error says that the SAM file does not have a header - check the file yourself - does it have a header?

BTW the `samtools view` step is not needed anymore, all samtools actions will auto detect and work seamlessly on SAM or BAM files. I also think that samtools now detects input stream, so specifying the `-` should not be needed in say sorting. The commands chain together much more elegantly now.
Question: Hi friends, I need your help please. I tried to adopt a tutorial in **human** (GSE27447) to my case in **citrus sinensis**, I searched a lot in google but I could not find fit equal to some parts and when I did these steps by removing parts I had doubts about, the result was a rma.txt which could not be read in some case (with ????) I do not know what I should type instead of... my query is **GSE67376 (citrus sinensis)**

```r
source("http://bioconductor.org/biocLite.R")

biocLite("GEOquery")
biocLite("affy")
biocLite("gcrma")
biocLite("hugene10stv1cdf") # changed to biocLite("citruscdf")
biocLite("hugene10stv1probe") # changed to biocLite("citrusprobe")
biocLite("hugene10stprobeset.db") ????????????
biocLite("hugene10sttranscriptcluster.db") ???????

library(GEOquery)
library(affy)
library(gcrma)
library(hugene10stv1cdf) # changed to library("citruscdf")
library(hugene10stv1probe) # changed to library("citrusprobe")
library(hugene10stprobeset.db) ?????????
library(hugene10sttranscriptcluster.db) ?????

setwd("C:/Users/Man/Desktop/New folder (3)")

getGEOSuppFiles("GSE27447")

setwd("C:/Users/Man/Desktop/New folder (3)/GSE27447")

untar("GSE27447_RAW.tar", exdir="data")

cels = list.files("data/", pattern = "CEL")
sapply(paste("data", cels, sep="/"), gunzip)
cels = list.files("data/", pattern = "CEL")

setwd("C:/Users/Man/Desktop/New folder (3)/data")
raw.data=ReadAffy(verbose=TRUE, filenames=cels, cdfname="hugene10stv1") ??????

data.rma.norm=rma(raw.data)
rma=exprs(data.rma.norm)
rma=format(rma, digits=5)

ls("package:hugene10stprobeset.db") #Annotations at the exon probeset level ?????????
ls("package:hugene10sttranscriptcluster.db") ??????????????

probes=row.names(rma)
Symbols = unlist(mget(probes, hugene10sttranscriptclusterSYMBOL, ifnotfound=NA)) ???????????
Entrez_IDs = unlist(mget(probes, hugene10sttranscriptclusterENTREZID, ifnotfound=NA)) ????????

rma=cbind(probes,Symbols,Entrez_IDs,rma)
write.table(rma, file = "rma.txt", quote = FALSE, sep = "\t", row.names = FALSE, col.names = TRUE)
```

Answer: ```r
biocLite("hugene10stprobeset.db") ????????????
biocLite("hugene10sttranscriptcluster.db") ???????

library(hugene10stprobeset.db) ?????????
library(hugene10sttranscriptcluster.db) ?????
```

These lines do nothing. I don't think that there is a package called citrus.db.

```r
raw.data=ReadAffy(verbose=TRUE, filenames=cels, cdfname="hugene10stv1") ??????
raw.data=ReadAffy(verbose=TRUE, filenames=cels, cdfname="citruscdf")
```

```r
ls("package:hugene10stprobeset.db") #Annotations at the exon probeset level?????????
ls("package:hugene10sttranscriptcluster.db") ??????????????
```

do nothing. this is to view the content/ functions of the package.

```r
Symbols = unlist(mget(probes, hugene10sttranscriptclusterSYMBOL, ifnotfound=NA)) ???????????
```

This is to get the gene symbols. I don't know how to get for citrus

```r
Entrez_IDs = unlist(mget(probes, hugene10sttranscriptclusterENTREZID, ifnotfound=NA)) ????????
```

This is to get the Entrez Ids. I don't know how to get for citrus
Question: do u have any recommendations for easier learning of enrichr?

Answer: Use it through [gget enrichr][1] for super easy use!

[1]: https://github.com/pachterlab/gget
Question: Hi,

I've done a model design and I hope someone can help out with my understanding of it!

I have an experimental setup that looks something like this:

**3 Time points (0hrs, 6hrs, 12hrs)**

**3 Different Conditions (Treatments A, B and C)**

So that makes **9 different combinations** of time points and treatments, each is in triplicate. There are therefore **27 Samples**.

My design formula is : `~ Treatment + TimePoint + Treatment:TimePoint`

My current understanding is this will give me small p-values of Treatment-specific effects over time?

I wanted to further refine this and look at a specific treatment and how it differs between two time points. So I used the following line:

```
foo <- list(c("TreatmentA.TimePoint12hrs"), c("TreatmentA.TimePoint0hrs"))
resMFType <- results(dds, contrast=foo)
```

Is this correct?

Thanks

Answer: "My current understanding is this will give me small pValues of Treatment-specific effects over time?"

I wouldn't say that, rather you're fitting with a model that can measure treatment-specific effects over time (i.e., time:treatment interactions as well as time-specific and treatment-specific effects). Whether the resulting p-values are small or not depends on whether there are any significant effects in the dataset at hand (yes, I'm being rather nit-picky here :) ).

BTW, you can shorten your design to ~Treatment*TimePoint.

The contrast you mentioned looks correct for looking at changes between 12 and 0hrs within treatmentA.

Anyway, you seem to have the correct design and know how to form the contrasts, so you should be good to go!
Question: Hi, I want to calculate RPKM values. I did DESeq and i got the results. My next step is i have to identify RPKM values for each gene. My pipeline is Bowtie>tophat>HTSeq>DESeq .plesase guide me through right path.

Answer: You can compute rpkm from a `DESeqDataSet` in the way indicated in @igor's answer. However, RPKMs should only be used for downstream analysis and not for testing differential expression. This is explicitly mentioned in the documentation of [DEseq2][1].

>In order to test for differential expression, we operate on raw counts and use discrete distributions as described in the previous Section 1.4. However for other downstream analyses – e.g. for visualization or clustering – it might be useful to work with transformed versions of the count data.

Same is said in the documentation for the older package [DEseq][2]:

>The count values must be raw counts of sequencing reads. This is important for DESeq’s statistical model to hold, as only the actual counts allow assessing the measurement precision correctly. Hence, please do do not supply other quantities, such as (rounded) normalized counts, or counts of covered base pairs – this will only lead to nonsensical results.

Other relevant discussion:

[Answer by Gordon Smyth about the use of RPKMs with voom, edgeR and DEseq2][3].

[1]: http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html
[2]: http://www.bioconductor.org/packages/release/bioc/html/DESeq.html
[3]: https://support.bioconductor.org/p/56275/#56299
Question: Hi,

I have a previously created Boxplots using normalized log2 data using `ggplot2` library in R. Now, I am trying to create the Boxplot using fold change values (example data below) in ggplot2. Could you please assist me how to do the same or provide any useful resources.

dput(head(Data))
structure(list(Genes = c("Gene A", "Gene B", "Gene C", "Gene D",
"Gene E", "Gene F"), Healthy_Control_p_val = c(1L, 1L, 1L, 1L,
1L, 1L), Healthy_Control_FC = c(1L, 1L, 1L, 1L, 1L, 1L), Condition_1_p_val = c(0.002651515,
2.48e-05, 0.73640094, 0.009591015, 0.501723078, 0.226255478),
Condition_1_FC = c(-1.114126532, 1.228203535, -1.054172268,
-1.255231765, 1.028948329, -1.086890065), Condition_2_p_val = c(0.043930415,
0.163740086, 0.798708578, 0.575282124, 0.292885948, 0.008037055
), Condition_2_Relapse_FC = c(-1.133444608, 1.163213331,
-1.195016491, -1.070544661, 1.051757975, -1.190927361), Condition_3_p_val = c(0.472872641,
0.62361238, 0.31757426, 0.481332218, 0.169341523, 0.781283534
), Condition_3_FC = c(-1.04283358, 1.04223539, 1.01103926,
-1.112300841, 1.074033832, -1.041999594)), row.names = c(NA,
6L), class = "data.frame")

Thank you,

Toufiq

Answer: You first need to Convert your data from wide to long format.

library("tidyverse")

df <- df %>%
select(!ends_with("p_val")) %>%
pivot_longer(ends_with("FC"), names_to="condition", values_to="fc")

> head(df, 5)
# A tibble: 5 x 3
Genes condition fc

1 Gene A Healthy_Control_FC 1
2 Gene A Condition_1_FC -1.11
3 Gene A Condition_2_Relapse_FC -1.13
4 Gene A Condition_3_FC -1.04
5 Gene B Healthy_Control_FC 1

Now that your data is in long format, you can make any sort of box plot you want. Here's an example.

ggplot(df, aes(x=condition, y=fc)) +
geom_boxplot() +
theme(axis.text.x=element_text(angle=45, hjust=1))

Plot of example data.

![enter image description here][1]

[1]: https://i.imgur.com/KYV7VRV.png
Question: Dear all,

I tried to open the example RNA-seq raw reads by using the simple command like head, but I got the error as shown below

$ head ~/RNA-seq/Liu2015/SRR1272191_IDCLV_1.fastq.gz
??˲?Ʋ%8?_Qf2K1?)?/? gp+̪~??{?֣??p_ `?p ?K???&?[?΀???Z??????|;?????1??????_?Yu?fU???????խ:???????T??????????????{x????????[?X?t?\?????xV;?c???@?{?tG?{{>??s?<$??mz??
}
??l???n?j??ț?\????+=??Z?eiSM??ö?-6v??????_??O?=lw???C?????#??lS??])??`?q?????x?a????l???n?????Z۫2o??}0??ME?Sx\?J/l?oYi??~????>؈ܠ??~??g钦?x??]?G??$w?o?????&{b?-Ά??T?eȮ?CW?NnV????T_?!???m????5?6?n??VU??t?w:<~?.??????Wz)r^i??Qta?=?M:?,_?w"?T?N;?=???ֶ?)L?z?Q?:???|?????e??????W?O????-v$.^?\?ˮ?ƕ(@?x?[???;S???	?+᳊???ޕ¦??`??#?=qZ??o?'z??`???"?e?m????}i?2?K?{Y6??X??oL???G?~?a???m?Oq?x\?|R?Ů?q?S??,~???*펝?x????>znaӽ+???
ۙ?o???*Gp

I didn't find any problem when I used this read to run FastQC and TrimGalore.

Could you please give me a suggestion, What should I solve this problem?

Thank you so much,

Best,

Kamoltip

Answer: Dear all,

Thank you so much for suggestions from everyone, I tried them all and I would like to inform everyone about my result.

For `$ zcat test.fastq.gz | head` it does not function in my case. It always shows the error. `zcat: can't stat: test.fastq.gz (test.fastq.gz.Z): No such file or directory`.

Anyway, `$ zless test.fastq.gz` and `$ zmore test.fastq.gz` work very well.

Thanks again for your kind help :)
Question: Dear all,

I am working on a method which combines the results of two procedures (limma, GSEA) to identify differentially expressed genes/pathways to obtain even better results using a simple machine learning approach.
For this, I need some kind of gold standard to compare my results with. So my question is: Does somebody know a dataset for which the results (up or down regulated genes/pathways) are known/sufficiently proved (e.g data set used by a paper)?
The data should be publicly available and microarray based.

Best regards,
Jan-Niklas

Answer: Hi Jan-Niklas,

please have a look at the SEQC experiments ([see here][1]) and the ABRF study (see publication [here][2]).

Both studies are based on the FDA-samples (A,B,C, & D) which are standardised mRNA-samples and have also been used in the MAQC studies for Microarray platforms. Both studies are publicly available at GEO.

Cheers,

Michael

[1]: http://www.fda.gov/ScienceResearch/BioinformaticsTools/MicroarrayQualityControlProject/#MAQC-IIIalsoknownasSEQC
[2]: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4167418/
Question: Hello everyone I'm looking for a function in python or Biopython that can calculate the `tetranucleotide frequency` of a given regions of scaffold.

The idea is that I have several **regions** and I want to identify possible `changes in nucleotide composition` that correspond to the an `endogenization` regions within my genome, for that I need to calculate theTNFs across regions for these contigs. I then need to calculate the Pearson correlation of these frequencies compared to the TNF of a set of the largest contigs in these genome assemblies (these contigs being probably really from the genome and were not endogenized).

Does someone know a such package in python?

Thanks you

Answer: [**CheckM**][1] can do what you need - see [**here**][2]:

checkm tetra seqs.fna tetra.tsv

You can use the frequencies to separate the sequences in a 2D plot using various dimensionality reductions methods such as PCA, tSNE (shown below) or UMAP.

![enter image description here][3]

[1]: https://github.com/Ecogenomics/CheckM
[2]: https://github.com/Ecogenomics/CheckM/wiki/Utility-Commands
[3]: https://i.ibb.co/xG0Zfmt/t-SNE-cluster-plot.png
Question: Before I report this as a bug can someone show me the proper way to add metadata using the `AddMetaData` function.

I first read a file that is new line delimited with the meta data labels I want using read table. The ftimecell_columns.txt file looks like this:

"12wks"
"12wks"
"12wks"
"12wks"
"12wks"
...

So I do this command to make the dataframe:

timecell_col <- read.table("~/Projects/FetalPancreas/timecell_columns.txt", header=FALSE, sep="\n")

and then this to AddMetaData:

AddMetaData(object = scfp, metadata = timecell_col, col.name = 'time.cell')

However, when I do:

head(x = scfp@meta.data)

All I see is the original meta data values:

orig.ident nCount_RNA nFeature_RNA
GSM2978830 scfetpan 382565 5080
GSM2978831 scfetpan 634726 4932
GSM2978832 scfetpan 565912 4501
GSM2978833 scfetpan 717152 5117
GSM2978834 scfetpan 207659 3508
GSM2978835 scfetpan 869148 5055

Am I reading the data in correctly?

Any help would be appreciated!

Very Respectfully,
Pratik

Answer: So figured it out. Months later, but nonetheless got it.

For reference: My column names were the GSM sample names because I downloaded the data from GEO, where each GSM####### corresponds to a cell (see image below:

![enter image description here][1]

Your gene expression matrix/data frame should be cells as column names, and genes as row names.

It may be obvious to some, but it was not obvious to me then. In order to add meta data. Your **meta data, data frame should have your row names as cells** and the **corresponding meta data information should be the row values, so it should look something like this**, and lastly **your column names should be what you want that meta data to be named** see below:

![enter image description here][2]

If your meta data, data frame is arranged in the way. You could simply do a modification of the default Seurat command like this by adding the `meta.data = yourmetadata` like so:

scfp <- CreateSeuratObject(counts = df, project = 'scfp', meta.data = SCFPmetadata, min.cells = 3, min.features = 200)

or the following after your make your Seurat object:

scfp <- AddMetaData(object = scfp, metadata = SCFPmetadata)

[1]: /media/images/7f0a6f5c-47f6-4b38-9267-61b235ee
[2]: /media/images/9ac4db2f-b63a-44b9-a85f-c8c127c8

Hope this helps!
Question: Hi all! I am having problems and I hope I can get some help from you.

I will explain my situation: I'm trying to perform a PCA analysis to see how different several bam files are. I'm using the next pipeline:

1.
Getting the accession files. I am using the R library "SRAdb", so I am getting 4 files in .sra format.
2. I use SRA-tools in order to convert the .sra file into .bam format
with the following code:

`sam-dump -r --min-mapq 25 $file | samtools view -bS > $file.bam`

3. Sort `samtools sort $file -o $file_sorted`

4. Index `samtools index $file_sorted $file_sorted.bai `

5. Compute a matrix to generate the PCA plot

` multiBamSummary bins -b $files.bam -o my/out/path --smartLabels -bs 10000 -p 2`

At this point I'm getting the following error:

The file < myfile > does not have BAM or CRAM format.

I haven't been able to trace the error, as any of the earlier steps reported any source of error. Any suggestions? (ideally I would like to skip the alignment step, I want to keep the file as original as possible)

- sra-tools --version 2.9.1_1
- samtools --version 1.9
- deeptools --version 3.3.0

Thanks before hand!!

Answer: Let us use one of the example accession numbers above (SRR2089860). These are single-end reads.

Your options are:

Use `fastq-dump` to dump the reads out in fastq format (remove `-X 5` for full set)

$ fastq-dump -X 5 SRR2089860
Read 5 spots for SRR2089860
Written 5 spots for SRR2089860

Use `sam-dump` to create fastq format files

$ sam-dump --fastq SRR2089860 | head -16
@HWI-D00473:169:HFK7WADXX:1:1101:1202:2011/1 unaligned
NGAGTCTATACTCGTTACATTCGCGTAACTCATTGTTAATCGCGAAGTTGA
+
#1=DDDDFGHHGHJJJIJIIJGIIJJIJHICGIIIJJJIJGIJEHJIGIIG
@HWI-D00473:169:HFK7WADXX:1:1101:1195:2074/1 unaligned
CTCGAACTCCTCGTAGTGGCGATTGTCGGTGCTGCCCACCAGGTCCACTGT
+
CCCFFFFFHGHHHJIJIJJJJHIIGGGHIECEHFHGIEFIGGJGHJIIGIG
@HWI-D00473:169:HFK7WADXX:1:1101:1230:2087/1 unaligned
TGCCGGGAATTGTACAGTGCTCAGCTTTATAGGACATTTCCAAACAGTTAT
+
BBBFFFF8FHHHHJJJIJJJJIJGJJIJFGJIFGIIIJJJIGIEIIIIJGG
@HWI-D00473:169:HFK7WADXX:1:1101:1222:2168/1 unaligned
CCGAGACTTGCCTGCTCACCAGCGAAGAGGGCGAGGAGCGTTTGACGGCCG
+
@@CDDADDHFHHHIIIIIHGGE
Question: Good night, could you please help me with the following question. Is it better to use logFC or logCPM to analyze my RNA-seq data between different treatments?
which is better to use p-value or FDR?
thanks for your help

Answer: Since logFC reflects the *difference* between your conditions, and that's what you're interested in, that is what you should pay attention to, and what will be meaningful to think about. LogCPM means very little in RNA Seq experiments. If you have RNA-Seq data, you're typically measuring thousands of genes, which means you're testing thousands of hypotheses, and so it is better to use FDR rather than simple p-values. Remember that for a p-value cut off of 0.05, you're essentially saying you are rejecting the null hypothesis (no difference between your conditions) in favor of the alternate hypothesis (there is an effect, a measurable difference between your conditions), and that if you're wrong (the null hypothesis is actually true), you would only see an effect as large as the one you measured 1 in 20 times. You can apply this same logic to all the genes you're measuring. Under the null hypothesis and a p-value cutoff of 0.05 you would expect a false positive 1 in 20 times, and since you're measuring thousands of genes, you can expect to have many genes pass a p-value threshold by chance simply because you are performing so many measurements. Thus you must "adjust" your p-values to account for this (calculate an FDR), which usually means inflating the p-values in some way. Genes with very low p-values will survive the adjustment (a very tiny number can still be very tiny even if multiplied by another number). (if you search here for FDR or adjusted p-value, you'll find better explanations by people who actually know what they're taking about). Also, make sure you read the [edgeR userguide][1], it's a good reference and explains what CPM is and isn't good for).

[1]: https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf