cannin/biostars_qa huggingface dataset:
Question: Hey guys:
I am doing **RNA-seq** analysis and it seems that the quality of my reads is not desirable.
Below is **a typical fastqc report** for my data.
I have read many tutorial about fastqc, from my understanding, **it seems that the 1-10 bp are adaptor sequences**. But in the adaptor content section, there is no waining.

----------

- I am wandering if my understanding is right?
- Should I use **trimmomatic** to cut adaptor sequences?

![Failed per base sequence content][1]

[1]: https://lh3.googleusercontent.com/Xp-_0jmKsw8p2MVrOx_X-hJbHIiKBfUFZ700lNT5nv0MKlma6DSgYonxLAbkHzkYNEaWHszNcvI6RjCPhk6YghTHfI7dT1_2GUyzPr4I5MdEnee1_qTSB46nbvwTJFh_TEASRgS1MEGIawbdc7zZpnWWd_bFL_CeNB0-MBk6PePyflSrMjeLZFlFu-1onTqcypuSMKeYw5Wty0c-pgPSCHWIHkTxrg3d0Y-XUhztxFRETL2sFQQXjtSGmZ5PhguduitzrZcnOijNlVK5KK_0hZ5Aqd-3pf7OHjJbRXfP02dxda5HQOrALIo-_DYoHszlbMDSFeLwSQerz80bPd2CNEsLmE79KsYBSqjOIeg_aH8npy70IlZc2pT34_pZul8pZapn_Yb_QN2tiazcaPPeYJvC-uVttpkoepsoaXT7Oz5lbB1OBkCrx1Ye3dyzpAx_bGhvdnnIXhOFYKbI44Nr1mh1ns7up0uj6uukZDqFOvRAkEFRNeNPHuqZlxGFfnCY1qAuM3Kqy3c0UjhVdDmFlTNVnR6wICnRR8v4EL7HxCbdPAnZ7JlHHzpaPHtKTxlMG80d51jSwWVLPa1GnaveRBsjfp3Ek5NVZ8bhv5UxhQXhPKRi8L3XBZ6-bg0t3MGomxWRSaE54EN6HD2GV1NKXE8muzTGcw=w1407-h809-no

Answer: The adapter starts at the 3' end of your reads, not the 5' (unless its an adapter dimer - i.e. no insert).

This is the result of random priming in RNA-Seq. I think [Biases in Illumina transcriptome sequencing caused by random hexamer priming][1] is the first paper on this. This represents real sequences. After alignment you can check the error rate in the reads, its only marginally higher 5' than in the rest of the read. And the fastqc report of the aligned sequences should show the same pattern.

[1]: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2896536/
Question: Hi All,
I have a question regarding trimming of adapters in NGS data. I have previously analyzed RNAseq data using brtim (https://www.sciencedirect.com/science/article/pii/S0888754311001339) without using adapter sequence. I am now analyzing smallRNAseq data, and I am using adapter sequence as adapter.fa using BBduk.sh from BBmap tools. Could you please clarify under what circumstances one would need to provide/know adapter sequence and when it is not necessary?

Answer: One case is where there is a kit/prep specific adapter that is being used. This may require special handling of the downstream data based on the instructions included in the kit.

If you had paired end reads with enough pairs having short inserts then you can detect them by doing this: `bbmerge.sh in1=r1.fq in2=r2.fq outa=adapters.fa`
Question: How to embed known cell information into Seurat object? Using AddMetaData?

Is this possible?

I read this documentation: https://www.rdocumentation.org/packages/Seurat/versions/3.1.4/topics/AddMetaData

However I need some clarification on how to go about doing this, please?

I have cell labels that are all aligned with the expression matrix such as this:

"12wks_fetal pancreas cell_acinar"	"19wks_fetal pancreas cell_beta"	"12wks_fetal pancreas cell_ductal"	"12wks_fetal pancreas cell_beta"	"12wks_fetal pancreas cell_acinar"	"22wks_fetal pancreas cell_acinar"	"12wks_fetal pancreas cell_acinar"	"12wks_fetal pancreas cell_alpha"	"12wks_fetal pancreas cell_ductal"	"12wks_fetal pancreas cell_acinar"	"14wks_fetal pancreas cell_endocrine.progenitor..."

I used these cell labels as cell names (column names) for the expression matrix. However that's not useful for seeing which cell is from what time point 12wks, 19wks, 14 wks, etc. on the UMAP for example.

Would there be a way to label the cells using the AddMetaData function so that when I visualize the data using UMAP, they are color coded based on the time of collection (ie. 12wks, 19wks, 22wks, etc.) or would there be a way I could also filter by both time of collection and/or cell type?

Is this possible in Seurat? How would I approach this?

I would really appreciate anyone's help.

Very Respectfully,
Pratik

Answer: Hope I got your question correctly, you can do everything to cell information via `Ident` function. Follow the manual to change the `Idents` of your cells. By this, you can replace the `active.ident` with what you want.
Question: **Hi everyone**

I would like to take a second opinion regarding a variant with DP =17. The variant is in a very high GC content area so the low coverage is reasonable.

I do hard filtering because I have few samples for VQSR, I thinks it is okay according to GATK hard filtering, the values are as follows:

AC=1;AF=0.500;AN=2;BaseQRankSum=0.892;ClippingRankSum=0.000;DP=17;ExcessHet=3.0103;FS=10.843;
MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=2.81;ReadPosRankSum=-1.534;
SOR=2.795 GT:AD:DP:GQ:PL 0/1:14,3:17:76:76,0,506

I am worried because **AD is 14,3** which means only 3 reads only support the ALT? is it safe to accept it as a het variant. The variant is very interesting based on our data, so I don't wanna discard without a strong reason.

**Any comments, how do you think? with many many thanks**

Answer: I cannot provide a detailed answer given the information but a few comments.

There are 2 most likely models:
1) The site is homozygous reference

Therefore, the 3 reads are potentially a) mismapped, b) weird duplicates or c) sequencing errors.

a) Did you use a mappability filter? i recommend Heng Li's method lh3lh3.users.sourceforge.net/snpable.shtml
b) Did you use rmdup?
c) if both a) and b) are satisfied, then they are sequencing errors, assuming a base quality of 38, the probability of 3 seq. errors occurring independently is:

(10^((-1*38)/10))^3 = 3.981072e-12

so one chance in 250 billion.

2) The site is heterozygous reference

It is possible that this is a genuine het site and you haven't sampled the other base, this can happen with probability:

dbinom(3,17,prob=0.5) = 0.005187988

so 1 chance in ~200 (assuming no mismappings or errors).

This is reflected in your PL field where the homo ref model has a score of 76, the het model a score of 0 (most likely) and 506 for the homo alt. While the hetero is more likely that the homo ref, it is not astronomically more likely. However, the homo alt is so unlikely that it can be safely discarded.
Question: Hi,

I've recently started using Ensembl for help in designing a gene panel for NGS.

Each time I select a gene on Ensembl it comes up with different available transcripts for that gene.

What defines a transcript? Is it literally just a different version of the gene in different individuals (but if so, how come some of the transcripts appear to be so different in terms of bp length and amino acid number)

Any help in explaining would be most gratefully received!

Answer: This is more of a comment, possibly a rant. In my opinion the concept of gene is obsolete and we would be better off if we ditched the concept of "gene" altogether.

Genes made sense when it was thought that there were these discrete units (genes) which produced each a single transcript and a single protein. However, a single gene can produce multiple transcripts and these can be very different one from the other (see for example the table of transcripts of the [Actin gene][1], there are several coding and non coding transcripts). So when we say gene X has a mutation it is not clear what we are referring to. How many transcripts are affected by this mutation? Are they coding or pseudogenes? In my opinion it would be simpler to think in "transcript space" and forget about genes.

Another way if seeing this is to consider that while transcripts exist, genes don't exist. When you do mRNA extraction you isolate transcripts, not genes, you can't isolate genes. You could have a restriction enzyme that cuts left and right of a "gene", but this is also not true. You have an enzyme that cuts at positions A and B, which happen to include a bunch of transcripts. If the definition of that gene (and transcripts) change, the enzyme still cuts there because it doesn't "see" genes, it sees DNA sequence.

I think "genes" hang around because they make some statements simpler ("Mutation at position A hits gene X" as opposed to "Mutation at position A hits an intron of transcripts X, an exon of transcript Y, and a UTR of Z"), but they are at best incomplete statements.

Any thoughts?

[1]: http://www.ensembl.org/Homo_sapiens/Gene/Splice?db=core;g=ENSG00000075624;r=7:5527151-5563784
Question: Hi All,

I have MiSeq PE250 reads for a viral vector sample. I trimmed the raw reads with Trim Galore (length >= 200 and Q >= 30). I aligned these reads using both `BWA MEM` and `Bowtie2 --local` and found a `AAA -> TTT` variant in both alignment files at different frequencies `9%` and `71%` respectively. I further analysed reads containing variant and found that the reads with variant had base quality call score less than 30 for TTT bases.

To check the effect of trimming on variant call and coverage, I further trimmed the trim galore trimmed reads to `remove all reads which had even a single base call quality less than 30`. I repaired these reads using `bbmap repair.sh` and aligned again using bowtie2. For this alignment the above variant was found at `<10%` variant frequency. This was consistent with variant frequency reported with BWA-MEM. The coverage was also affected with almost `6.7%` of the bases with `coverage less than 10X` for the super trimmed reads vs `0.1%` for trim galore trimmed only reads.

I used freebayes for variant calling with same parameters (Min Coverage 10, Min Alternate Fraction 0.01, Min Alternate Count 4) for both BWA and Bowtie2 aligned reads.

1. Why was there such a large difference for variant frequency for Bowtie2 and BWA MEM?
2. The coverage difference between trim galore trimmed reads and super trimmed reads is vast. Should I still use super trimmed reads (Trim Galore + Further trimming) as it gives correct variant frequency.

The fastqc result for the trim galore trimmed reads and super trimmed reads were as follows:

**Trim Galore Trimmed Reads**

**Super Trimmed Reads**

The fastqc per base faultily fails for the super trimmed reads indicating a wrong Illumian Phred Score encoding version. Why did this happen?

Thanks!

Answer: >The fastqc per base faultily fails for the super trimmed reads indicating a wrong Illumian Phred Score encoding version. Why did this happen?

Because with this line:

```
awk '{unit=unit $0 ORS} NR%4==0{if (/^[?@ABCDEFGHIJK]+$/) printf "%s", unit; unit=""}' trimmed.R1.fq > super_trimmed.R1.fq
```

you eliminated any reads containing any bases with a score below `?`, which leads `FastQC` to the assumption that you're looking at files where the Phred scores have been encoded the Illumina1.3-way

![](https://i.ibb.co/nPYvV7Q/Screen-Shot-2019-04-01-at-4-56-17-PM.png)
Question: I have an ExpressionSet, combArr.eset, that contains Affymetrix Human Genome U133 Plus 2.0 Array gene expression data. I also have a matrix, exprs_combArr.eset, that contains the expression data for combArr.eset. The row names for exprs_combArr.eset are the hgu133plus2 probe IDs, and the column names are the Gene Expression Omnibus sample IDs.

I'd like to change the probe IDs in exprs_combArr.eset to their corresponding gene symbols, but I'm not sure how to do this.

I've retrieved the hgu133plus2SYMBOL object from the hgu122plus2.db package, which contains the mappings between manufacturer identifiers and gene abbreviations, but I'm not sure how to use this to change the probe IDs in exprs_combArr.eset.

I'm new to R and Bioconductor, and programming in general, so any help would be appreciated.
Thank you!

Answer: Haha, this came up in a biofinormatics clinic I did on Friday. I've written a rather terse version of how to add gene-level annotations to an Affymetrix ExpressionSet [here][1]. You basically have to find an appropriate annotation package (probably hgu133plus2.db) and add that info to the featureData slot of your ExpressionSet (don't monkey with the probe id-rownames)

[1]: http://biolearnr.blogspot.co.uk/2017/05/bfx-clinic-getting-up-to-date.html
Question: Hello,

I know the following question have been asked many times by now, but the suggested solutions were not helpful.
I have a set of taxIDs and I want to convert them to taxonomies, but I do not want to download any extra databases. Is there any simple way to do this via terminal?

Answer: Another solution using [Entrez Direct][1]

efetch -db taxonomy -id 387462,2594474 -format xml | \
xtract -pattern Taxon -first TaxId -element Taxon -block "*/Taxon" \
-unless Rank -equals "no rank" -tab "," -sep "_" -element Rank,ScientificName

387462 superkingdom_Eukaryota,kingdom_Metazoa,phylum_Mollusca,class_Gastropoda,subclass_Caenogastropoda,order_Littorinimorpha,superfamily_Stromboidea,family_Strombidae,genus_Strombus
2594474 superkingdom_Bacteria,phylum_Firmicutes,class_Clostridia,order_Clostridiales,family_Peptococcaceae,genus_Dehalobacter

[1]: https://www.ncbi.nlm.nih.gov/books/NBK179288/
Question: Hi,

I encountered a strange issue while reading in a data table from `txt` format. If I read it from `txt` by `read.table` it does not include all rows but if I convert to `csv` and read it with `read.csv` its perfect. Does someone know the issue or is it my code?

[Here][1] is the file.

```r
test <- read.table("./Annotations/all_genes_pombase.txt",
+ header=T,
+ sep="\t",
+ row.names=1,
+ stringsAsFactors = F)
> dim(test)
[1] 4533 7
> str(test)
'data.frame':	4533 obs. of 7 variables:
$ name : chr "SPAC1002.01" "pom34" "gls2" "taf11" ...
$ chromosome : chr "I" "I" "I" "I" ...
$ description : chr "conserved fungal protein " "nucleoporin Pom34 " "glucosidase II alpha subunit Gls2 " "transcription factor TFIID complex subunit Taf11 (predicted) " ...
$ feature_type: chr "protein_coding" "protein_coding" "protein_coding" "protein_coding" ...
$ strand : int 1 1 -1 -1 -1 -1 -1 -1 -1 -1 ...
$ start : int 1798347 1799061 1799915 1803624 1804548 1807270 1807996 1809480 1811408 1813740 ...
$ end : int 1799015 1800053 1803141 1804491 1806797 1807781 1809433 1811361 1813805 1815796 ...
```

[1]: https://drive.google.com/file/d/0Bxl0rZon0OuzR1VzN3cyMlhkUUk/view?usp=sharing

Answer: Use the argument `quote = ""` inside read.table.

read.table("your_file", quote="", other.arguments)

**Explanation**:
Your data has a single quote on 59th line (`( pyridoxamine 5'-phosphate oxidase (predicted)`). Then there is another single quote, which complements the single quote on line 59, is on line 137 `(5'-hydroxyl-kinase activity...)`. Everything within quote will be read as a single field of data, and quotes can include the newline character also. That's why you lose the lines in between. **`quote = ""` disables quoting altogether**.

There are other more instances where this 'quoting' happens again. One way to know how many fields read.table sees in every row is by using `count.fields`

num.fields = count.fields("all_genes_pombase.txt", sep="\t")

Now look at the variable `num.fields`, there will be a lot of NAs, the lines which are not read correctly by `read.table`

**The problem doesn't arise with read.csv** because the [quoting defaults][1] are different in read.table and read.csv, due to some reason really unknown to me!

read.table: quote = "\"'"
read.csv: quote = "\""

PS: The best way to avoid the reading file nuisance of *read.table* is to use *fread()* from data.table package. The side benefit is that it's blazing fast for large files and it guesses the field separator automatically. See my earlier post: https://www.biostars.org/p/221009/#221032

[1]: http://stat.ethz.ch/R-manual/R-devel/library/utils/html/read.table.html
Question: I find the way to get the identity of hsp ,but no alignment.One alignment can have multiple hsp, that means the score(identity) of hsp is no equal to alignment. when i do blast in webpage, i always get scores and identities of alignments. Is that any logical problem i made, who can explain this?

Answer: BLAST is a local alignment tool - it reports segments of the sequence, called HSPs, that align and produce the best scores. Each HSP in a match is a segment of query matching to a segment of the subject producing a high score. You can either use the meaningful maximum score or use a score totaled across all HSPs.

BLAST web UI has the same option - for identity percentage in case of multiple HSPs, it gives you the minimum identity% of the HSPs.
Question: Hi, everyone,

As my title describes, where to download *Arabidopsis* gene expression data in batch? I know some sites, for example, eFP browser, AtGenExpress Visualization Tool (AVT), but they only browser one gene each time. I would like to get global gene expression map. Could anyone help? Thank you very much!

Answer: Hello biolab,

If you are thinking about array-based gene expression, you could download in batch at [Araport][1]. There are at least two ways to do so:

1. If you don't have a specific experiment in mind, use this pre-defined template (Click [Here][2]). You could select the expression values of a list of genes. There are some demo lists available for you to get started. In order to use a costumed list you need to create a list by uploading the gene IDs (AGIs) [here][3]. The output of this query result is a interactive table that allows you to filter the experiment type, expression ration, etc.
2. If you are interested in a particular experiment condition, use this template instead (Click [here][4]). There is a drop-down menu allowing you to choose a specific experimental setting as well constrain the expression signal, ratio, and p-value. The result table would display all the genes met the criteria in that experiment.

Keep in mind that you would need to register and log in to permanently save the result. I hope this has addressed your question.

[1]: https://www.araport.org
[2]: https://apps.araport.org/thalemine/template.do?name=Gene_Expression&scope=all
[3]: https://apps.araport.org/thalemine/bag.do
[4]: https://apps.araport.org/thalemine/template.do?name=ExperimentCondition_gene&scope=all
Question: Hi,

I do have two groups of fastq files, I need to merge R1 reads of one file whose has same name at the beginning with the other R1 files reads.

Like Soil-13 is similar in two RI reads files,

I do have multiple pair-end read files like this that I need to merge into one. similarly, I need to do for R2 reads.

at the end I want to I have two soil-13 fastq files, one for R1 reads and other for R2 read. Like this, I need to do with multiple files.

Soil-13_S4_L001_R1_001.fastq

Soil-13_S4_L001_R2_001.fastq

Soil-15_S5_L001_R1_001.fastq

Soil-15_S5_L001_R2_001.fastq

Soil-13_S62_L001_R1_001.fastq

Soil-13_S62_L001_R2_001.fastq

Soil-15_S72_L001_R1_001.fastq

Soil-15_S72_L001_R2_001.fastq

Kind Regards

Answer: You can just concatenate them with cat:

cat Soil-13*_R1_*.fastq > Soil-13_R1_001.fastq
cat Soil-13*_R2_*.fastq > Soil-13_R2_001.fastq
etc.

You can generate a for loop to run through all libraries (bash):

for l in $(ls *_R1_*.fastq | cut -d "_" -f 1 |sort |uniq); do cat ${l}*_R1_*.fastq > ${l}_R1_001.fastq && cat ${l}*_R2_*.fastq; done
Question: Assume we identify - by RNA-seq, tiling arrays, by prediction - possible candidate regions for non-coding, small RNAs. I wish to verify and predict the function of as many RNAs as possible by computation before going to the lab. One could use eg. Rfam to find similar sequences, after that we are left with more than 90% that have no match. One could predict the 2D structure using eg. [RNAfold][1], compare that using [RNAforester][2]. But that does not get me even close to a function prediction. Do you have experience with other tools or a better computational pipeline that gets more information out of the ncRNA candidates, possibly even something specific to bacteria.

[1]: http://www.tbi.univie.ac.at/~ivo/RNA/RNAfold.html
[2]: http://bibiserv.techfak.uni-bielefeld.de/rnaforester/

Answer: How are you doing your comparison to Rfam? If you're just running rfam_scan.pl or the CMs then this may not give you the results you really want. See the recent [paper by Kolbe & Eddy][1] to hear more about the limitations of CMs on truncated sequences. Could explain some of your lack of sensitivity vs Rfam. However there is a terrifying number of ncRNAs not yet covered by Rfam.

I'm not sure clustering RNAfold predictions with RNAforester will tell you much either. [Locarna][2] and [CMfinder][3] have been used previously to cluster many ncRNA predictions.

[1]: http://www.ncbi.nlm.nih.gov/pubmed/19304875
[2]: http://www.ncbi.nlm.nih.gov/pubmed/20444875
[3]: http://www.ncbi.nlm.nih.gov/pubmed/16357030
Question: Hi All,

I have a table with six columns (in the following).

![enter image description here][1]

[1]: http://uupload.ir/files/xm5g_taible.png

So, I want to delete SNPs that are less than number 10 in the SNP column.

What is the best idea?

Best Regard

Mostafa

Answer: awk !

a cmdline such as below one should do the trick:

awk ' $4 >= 10' > new_file

and if you want to (additionally) apply calculations on a certain column:

awk ' $4 >= 10' | awk '$2=$2+2000' > new_file
Question: I am writing a script to get SNP positions from dbSNP by using NCBI E-Utils. When I use esearch and efetch to retrieve the positions, I can only access positions from the newest genome build. Since I want to use GRCh37, this causes problem. I was wondering if you knew of a way of asking for a specific assembly in E-Utils.

Answer: dbSNP does not support searches for older assembly coordinates.

Depending on exactly what you are trying to do, you can probably use lift-over or remap to first convert your coordinates.
Question: I simulated reads for a list of variant sites using GATK SimulateReadsForVariants tool. From that, I get a bam file output. Next I create a pileup using samtools mpileup.

samtools mpileup -t DP,DPR,DV,INFO/DPR -vuf GRCh37.fa -l snp.file.vcf simreads.bam > simreads.raw.vcf

The problem is that the output is not giving me correct counts at the variant sites. Here's an example of the first two SNP sites from the pileup.

```
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA12878
1 837214 . G 0 . DP=20;I16=0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;QS=0,0;MQ0F=0;DPR=0,0 PL:DP:DV:DPR 0,0,0:0:0:0,0
1 851390 . G 0 . DP=20;I16=0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;QS=0,0;MQ0F=0;DPR=0,0 PL:DP:DV:DPR 0,0,0:0:0:0,0
```

So why would all the counts (with the exception of the DP tag) be 0? I opened up the bam file in IGV to check the first couple of sites and the reads do pileup there as they should.

Answer: I figured out the solution myself. I had simulated more error in my reads than I initially thought. For anyone else experiencing a similar issue, check the `-Q`, `--min-BQ` mpileup flag. Samtools automatically does not count reads in the other tags unless they are "high quality" (i.e. have a minimum base quality of 13). Depending on your experiment, if you're only interest is getting all the reads counted then consider using `--min-BQ 1`.
Question: Hi,

I have some RNA-seq samples that I want to normalize and then output RPKM expression, but I am unsure how to do this.

This is my pipeline so far:

Normalise raw read counts with TMM in edgeR

expr <- DGEList(counts=data, group=conditions)
expr <- calcNormFactors(expr)

output:

$samples
group lib.size norm.factors
Sample1 F 19770521 1.0462660
Sample2 F 17970679 0.8794805
Sample3 F 19184265 1.0573665

QUESTION: How do I get normalized raw read counts from this? Do I multiply the read counts by the norm.factors?

QUESTION: Ultimately, I want to end up with RPKM values for each gene in each sample. I know I can use the rpkm() function below in edgeR

expr_norm <- rpkm(expr, log=FALSE,gene.length=vector)

but is expr the output from calcNormFactors or something else?

Thanks for your help!

A

Answer: The library size normalized counts are made by dividing the counts by the normalization factor (you'll note that the larger libraries have larger normalization factors, so if you multiplied things you'd just inflate the difference in sequencing depth).

For the rpkms, just do `rpkm(expr, gene.length=vector)`, since it can take your `DGEList`, (this is explained in `help(rpkm)`).