cannin/biostars_qa huggingface dataset:
Question: There are many experiments where samples are stored in 384-well plates (or 96-well plates etc.) and a measurement is made on each sample. Typical technical biases may affect spatial regions of the plate, for instance the borders (evaporation, ...), or large bands when the thermocycling (in PCR) is preformed by separate independent controllers, etc. In some more complicated experiments, each well may contain a NGS library, for which multiple quality control scores can be computed.

I searched for `R` functions for pretty printing of values arranged in the geometry of a multiwell plate and I was surprised that I did not find anything obvious. Perhaps because the solution is too trivial ? Still, I would be happy to use a well-thought packaged function instead of a quick-and-dirty ad-hoc home-made wheel reinvention.

Here is an example below. The problems with this script are that 1) it requires flipping some levels that should stay as they are otherwise and 2) SVG or PDF version had a very ugly appearance as each square was strongly smoothed...

*My question is whether there are good packaged tools (in R) to do such plots and more operations on representations of multiwell plates.*

plate <- data.frame( Row = rep(LETTERS[1:16], 24)
, Col = unlist(lapply(1:24, rep, 16)))

plate$Row <- factor(plate$Row, levels=rev(levels(plate$Row)))

set.seed(1)
plate[001:096, "val"] <- rpois(96,10)
plate[097:384, "val"] <- rpois(96, 6)

png("plate.png", h=400,w=600)
ggplot2::ggplot(plate, ggplot2::aes(Col, Row, fill=val)) + ggplot2::geom_raster()
dev.off()

![384-well plate with random values][1]

[1]: https://gist.githubusercontent.com/charles-plessy/6325a878429fef8aae6d9ed25f0054ef/raw/f56cb303524f8c943d86a62583cadea17ea4959f/plate.png

Answer: [HTqPCR][1] has a function named plotCtArray and another named plotCtCard. They are not so stylish, but you can avoid having to create a new function.

[1]: http://www.bioconductor.org/packages/release/bioc/html/HTqPCR.html
Question: Is there a way to tell blast to output the unmapped reads into a separate file or at least the name of the unmapped sequences?

Answer: I assume you have used blast to align some reads "reads.fa" onto a reference "ref.fa", and have a blast report in "out.bls" ?

% formatdb -i ref.fa -p F
% blastall -p blastn -i reads.fa -d ./ref.fa -o out.bls

As RM answered, you can (hackily) get the read IDs that got no hits. I have added extra code to ensure a clean list of IDs only into "nohits.ids":

% grep -B5 "***** No hits" out.bls | grep '^Query=' | sed 's/^Query= //' > nohits.ids

Now if you want to get back the sequences listed in nohits.ids from reads.fa, you can use this trick:

% formatdb -i reads.fa -p F
% fastacmd -d ./reads.fa -i nohits.ids -D 1 -o nohits.fa

Good luck!
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: `seqkit seq -r seq.fa > rev.fa`.

Install `seqkit` from [**here**][1].

[1]: https://bioinf.shenwei.me/seqkit/download/
Question: Hi,

So I am trying to add the gene ID to my RNA-Seq data. So what I did was use Salmon to quantify my reads but instead of doing some DE analysis, I want to try and us an alternative program to try and do a little machine learning. So I really just need add the gene symbol next to the transcript ID so I can easily identify the gene instead of having to look up the transcript ID in ensembl

I know there are some packages out there that can help with this like `org.Cf.eg.db`. I can't seem to figure out how to make the package work with transcript IDs though. I am inexperienced with the package so that is most likely my issue. Here is an example of what my data looks like after it was quantified.

> head(test)
Name Length EffectiveLength TPM NumReads
1 ENSCAFT00000034820.1 957 736.829 1309.272 43423.000
2 ENSCAFT00000034824.1 1044 823.630 1001.516 37129.000
3 ENSCAFT00000034830.1 1545 1324.630 3796.436 226357.000
4 ENSCAFT00000034833.1 684 464.046 8086.686 168910.000
5 ENSCAFT00000034835.1 204 50.476 4033.303 9163.596
6 ENSCAFT00000034836.1 681 461.059 9748.035 202300.391

I really just want to add the gene symbol so the data looks something like this

> head(test)
Name Length EffectiveLength TPM NumReads Symbol
1 ENSCAFT00000034820.1 957 736.829 1309.272 43423.000 ABC
2 ENSCAFT00000034824.1 1044 823.630 1001.516 37129.000 ABD
3 ENSCAFT00000034830.1 1545 1324.630 3796.436 226357.000 ABE
4 ENSCAFT00000034833.1 684 464.046 8086.686 168910.000 ABF
5 ENSCAFT00000034835.1 204 50.476 4033.303 9163.596 ABG
6 ENSCAFT00000034836.1 681 461.059 9748.035 202300.391 ABH

Answer: you can use R package "biomaRt" to annotate you transcript id to gene name. see if the below code works

library( "biomaRt" )
mart = useMart('ensembl')
# list all the ensembl database of organisms
listDatasets(mart)
#choose database of your interest ; in this case its "cfamiliaris_gene_ensembl" I guess
ensembl = useMart( "ensembl", dataset = "cfamiliaris_gene_ensembl" )
# choose attributes of your interest
listAttributes(ensembl)
gene <- getBM( attributes = c("ensembl_transcript_id","external_gene_name"),values = test$Name,mart = ensembl)
#Macth your transcript id with ensembl_transcript_id
id <- match(test$Name , gene$ensembl_transcript_id)
#Add Gene symbol column in your data frame
test$Symbol <- gene$external_gene_name[id]
head(test)
Question: Hello,

I currently have two GTFs which are incomplete. What I mean by this?. For example, in one of them I have 5000 genes and in another I have 7000, however, there are 3000 genes that coincide (they are the same). Therefore, I would like to create a GTF with the 5000 + 7000 genes, but without duplicating the 3000 identical genes, so I would have in the end a GTF with 9000 genes (not 12000).

What could I do? I think both Cuffmerge and Cuffcompare are useless in these cases.

All the best

Answer: To give more precision to the @genomax comment:

`agat_sp_compare_annotations.pl` will tell you the differences.

Then you can choose to **complement** the annotation using `agat_sp_complement_annotations.pl`, that means you use one annotation as reference and you will pick from the second annotation only non-overlaping predictions and add them to the reference annotation (i.e. for overlapping loci only the gene prediction from the reference is kept)

Or to **merge** annotation using `agat_sp_merge_annotations.pl`, that mean you take non overlaping locus from both annotations and you merge the annotation when they are overlapping (overlapping genes will be merged as one locus and different mRNA will be seen as isoforms, duplicates mRNA will be removed to keep one)
Question: Hi
I am new in python and want to see how this can be done in python (I can do this in R). I have a text file `myfile.txt` with one column and thousands of rows as shown below. I want to convert this to fasta `result.fasta` format as shown below. How can I do this in python?

`myfile.txt`

ATGTGTGGTTTTCCCCC
ATTGGCGGGGTTTTTCAGGGG
ATGGGGGGGCCCCCCCCAAAAAA
TTGGTGGGGGGGGGGGGAA

`result.fasta`

>1
ATGTGTGGTTTTCCCCC
>2
ATTGGCGGGGTTTTTCAGGGG
>3
ATGGGGGGGCCCCCCCCAAAAAA
>4
TTGGTGGGGGGGGGGGGAA

Answer: Something a bit simpler and should work in Python 2 and 3:

#!/usr/bin/env python

import sys

c = 1
for l in sys.stdin:
sys.stdout.write(">%d\n%s\n" % (c, l))
c += 1

Usage:

$ convert.py < in.txt > out.fa
Question: I just got my RNA seq data.

CZ1_S6_L001_R1_001.fastq

CZ1_S6_L002_R1_001.fastq

CZ1_S6_L003_R1_001.fastq

CZ1_S6_L004_R1_001.fastq

I only know there are from 4 different lanes during sequencing.

What are they and do I need to combine these fastq files into one? And how?

Thanks~~~

Answer: It appears to be the same sample run in 4 separate lanes. You can `cat` the files together into one or process them independently (giving you a way to parallelize). You can merge the 4 sample bam files into one (and then sort) at the end.
Question: Dear all,

I want to do differential expression at **transcript level** using an RNA-seq dataset that consists of **paired samples** (patient 1 control/treatment, patient2 control/treatment, etc.).

I have not found methods other than Ballgown and DEXSeq that are able to handle paired designs. I have also tried the option of obtaining **transcript counts** and use them in edgeR. However I am concerned about the fact that edgeR is not particularly adapted to do analysis on transcript counts.

Does anyone has any advice on how to go on this problem? Or if there is a way to assess my results with the latter approach I described?

Many thanks,
Maria

Answer: Perhaps not the solution you were looking for, but software like kallisto and Salmon (and others) generate estimates of raw counts per transcript. You could use those estimates in edgeR or DESeq2.
Question: I am trying to get the mean depth coverage for bam files with bedtools.

So far I've been using the command:

```bedtools coverage -hist -abam bamfile.bam -b targets_exome.bed > bamfile.bam.hist.all.txt```

But I find myself with a result similar to:

1	9998	10030	HWI-H217:64:C48Y8ACXX:4:1110:12530:87823/1	1	-	9998	10030	0,0,0	1	32,	0,	0	32	32	1.0000000

1	9998	10030	HWI-H217:64:C48Y8ACXX:4:1110:12530:87823/2	1	-	9998	10030	0,0,0	1	32,	0,	0	32	32	1.0000000

1	9999	10096	HWI-H217:64:C48Y8ACXX:4:1303:11532:70370/1	0	-	9999	10096	0,0,0	1	97,	0,	0	97	97	1.0000000

1	9999	10096	HWI-H217:64:C48Y8ACXX:4:1303:11532:70370/2	0	-	9999	10096	0,0,0	1	97,	0,	0	97	97	1.0000000

where the different features are overlapping (see #1 and #3).

I want to compute average coverage, defined as the number of reads covering the captured coding sequence of a haploid reference. How to do it in this case?

Answer: To get the mean depth of coverage for each interval specified in your BED file, use:

bedtools coverage -a BED -b [BAM] -mean > MeanCoverageBED.bedgraph

If you wanted to get the overall mean, then just pipe this into awk and get the average of the average, something like:

bedtools coverage -a BED -b [BAM] -mean | awk '{total+=$4} END {print total/NR}'

There are other commands with bedtools that I use as standard in my pipelines:

# output depth of coverage for all regions in the BAM file

*NB - sequential positions at the same read depth are merged into a single region*

bedtools genomecov -ibam [BAM] -bga -split > CoverageTotal.bedgraph

# output the per base read depth for each region in the BED file

bedtools coverage -a [BED] -b [BAM] -d > PerBaseDepthBED.bedgraph

# get percent genome covered

zero=$(bedtools genomecov -ibam [BAM] -g [ref.fasta] -bga | awk '$4==0 {bpCountZero+=($3-$2)} {print bpCountZero}' | tail -1)
nonzero=$(bedtools genomecov -ibam [BAM] -g [ref.fasta] -bga | awk '$4>0 {bpCountNonZero+=($3-$2)} {print bpCountNonZero}' | tail -1)
percent=$(bc <<< "scale=6; ($nonzero / ($zero + $nonzero))*100")
Question: Dear Biostar users,

I have been working with metagenomics data.

In order to perform binning I need to create depth file.

- For this I first ran `bbwrap` with `spades` contigs as reference and samples as input reads.

```
for i in ${finallist}
do
bbwrap.sh \
ref=spades/contigs.fasta \
in=${i},${i%1*}u.final.clean.fq \
in2=${i%1*}2.merged.final.clean.fq,NULL \
t=24 \
out= ${i%1*}sam.gz \
covstats= ${i%1*}coverage
done
```

- After this I did `samtools view` with follow command:

```
samtools view \
-S \
-b \
-u \
${i%1*}sam.gz > \
${i%1*}bam
```

- Followed by sorting the generated bam file using `samtools sort` function with following command:

```
samtools sort \
-m 24G \
-@ 3 \
${i%1*}bam \
-o ${i%1*}sorted.bam
```

**During the sorting I getting strange character error messages or terminal with empty file in the end

I tried to solve the problem by checking the header of my reference and bam file and it looks the same.**

**Here are first five lines of my bam file:**

```
NS500633:12:H7JF2BGXX:1:11101:4115:8221 83 NODE_179_length_41253_cov_14.682023 5385 44 121= = 5235 -271GTGGAAAATGAAAAATATATTGATATCCTATGGGAAGGCTATGGCAAGGGGGAAACGGTGACGATTGGCGAAGATAAGGAATTCCTCATGGCAGATGCAAAGTCAAAGGATGGAATGGAAC FFFF

```
>NODE_1_length_199624_cov_13.770641
>NODE_2_length_193066_cov_15.912808
>NODE_3_length_183817_cov_14.803240
>NODE_4_length_163959_cov_11.179160
>NODE_5_length_128783_cov_13.807804
```

I am really unable to understand where is the problem coming from. It would really great it someone can help with this.

Thanks in advance!

Cheers,

AR

Answer: The issue is resolved. I had installed the samtools through conda and it did not install HTSLIB and that was the issue. Now I have installed samtools locally and it ran normally.

Thank you helping me think beyond the script.

Cheers,
AR
Question: I have assembled a bacterial genome using spades.
The coverage among the contigs is very heterogeneous, while some have more than 30 fold of coverage, other have less than 5 fold. Also, there are a lot of contigs with a length less than 500 bp.
I'm aware that I must filter contigs having low coverage and small length (they probably are spurious sequencing products), however what would be reasonable criteria for filtering them?
I have found scripts that remove contigs having less than 500 pb and 2 fold coverage. I believe this setting is too permissive, since most of the contigs of my assembly significantly deviate from it.

EDIT:
![Here is the distribution of log10 contig length:][1]

![Here is the distribution of contig coverage:][2]

As we can see, most of contigs have less than 1000 bp, Most of these contigs have less than 5 fold of coverage. My plan is to filter contigs having less than 1000 bp and 5 fold of coverage. After, I am going to check the completeness of the filtered genome sequence using checkm (contrasting to the unfiltered genome).

EDIT 2:
![Length vs Coverage][3]

[1]: https://image.ibb.co/iLswOb/semia4088_log_leng.png
[2]: https://image.ibb.co/i7sMOb/semia4088_cov.png
[3]: https://image.ibb.co/jVC6PG/semia4088_len_vs_coverage.png

Answer: If you're afraid of contamination or spurious sequencing products, I suggest to align by BLAST all of the contigs to some reference (or even to the NCBI NT or NCBI NR databases entirely), to decide which of the contigs are "Ok" (i.e. have matches to the species you are interested in). Then you can plot the distribution of coverages of good and bad contigs and define some threshold that separates them good enough
Question: Say I want to parse out all split reads. These are reads with a clipped position and a secondary mapping.

How do I identify these reads in samtools? Is there a flag to find reads that are split with a secondary mapping.

Once I have these reads, how I can I tell which chromosome they are mapped to?

Here's an example of what I'm looking for:

The blue read should have a clipped signal, I want to parse out these reads that have the alternate mapping position on the same chromosome.

Thanks!

Answer: samtools view -f 256 Input.bam | awk '$6 ~/S/ && $7 == "=" {print $0}' > Secondary_clipped.sam

Please debug and use it.

PS: post modified.
Question: Hi,

What are people using to search for raw fastq files in GEO or SRA made available from past publications for a specific cell?
I am currently using SRA search database with key words, specifying in the search term:

GM12878[All Fields] AND cluster_public[prop]

However this misses some downloadable data.

There is a nice post here (https://www.biostars.org/p/111040/) to retrieve these files using fastq-dump, but this assumes the GEO accession number, or the SRA project page, is already known.

Is there another searchable engine that helps to search for these more efficiently that people are currently using?
I see there is an R package called GEOquery to help with GEO searches - I haven't looked at it yet but if it is possible to limit by only certain files and cell types, this may be the best option.

Thanks much!

Answer: You can try `sra-explorer` ( https://sra-explorer.info/ ) from Phil Ewels.
Question: Hello

I wonder if there's a way to get the information provided by http://useast.ensembl.org/Homo_sapiens/Tools/IDMapper programmatically, through biomaRt or some other software.
Maybe this information is stored somewhere in ensembl ftp release?

Specifically, I have a list of human ensembl IDs for build 84 and for build 93, and I need to reconcile them.

Answer: The web tool is the front end to a [script using the perl API available in the Ensembl GitHub repository][1].

[1]: https://github.com/Ensembl/ensembl-tools/tree/release/99/scripts/id_history_converter
Question: I'm trying to use SingleR to annotate UMAP plots in Seurat. I've been following this very helpful tutorial here (http://bioconductor.org/books/devel/SingleRBook/pancreas-case-study.html) and I see that SingleR can be used to annotate clusters in Scran but I haven't been able to figure out how to do the same thing in Seurat. Would someone be able to help me do it? I've pasted a portion of my code below to give a sense of what I want to achieve and what I have been trying.

```
library(Seurat)
library(DropletUtils)
library(SingleR)
library(celldex)

# Write seurat_integrated to 10X file
write10xCounts(x=seurat_integrated@assays$RNA@counts, path='seurat_integrated_10
X/')

# Annotate with SingleR
new.data <- Read10X("seurat_integrated_10X")
ref.data <- MouseRNAseqData()

predictions <- SingleR(test=new.data, assay.type.test="integrated",
ref=ref.data, labels=ref.data$label.main,
de.method="wilcox")
table(predictions$labels)

# Annotate UMAP clusters
seurat_integrated$celltypes <- predictions$pruned.labels
DimPlot(dummy_seurat_integrated, reduction="umap", split.by="celltypes")
```
I've had trouble pasting it in here but in the link, there is a SingleR-labelled TSNE plot achieved using Scran. I'm trying to do the same thing, but using Seurat.

Answer: Provide the `clusters` parameter to your `SingleR` call and basically proceed as you've shown. I'd recommend reading the `SingleR` function documentation, i.e. `?SingleR` from within R.
Question: Hello,

I have a data in the form of a dataframe that I downloaded from [GTEx Portal][1]. It contains RNASeq gene read counts used in their study.

```
> dim(expr.df)
[1] 55993 2921
> expr.df[1:10,1:2]
GTEX-N7MS-0007-SM-2D7W1 GTEX-N7MS-0008-SM-4E3JI
ENSG00000223972 0 0
ENSG00000227232 158 166
ENSG00000243485 0 0
ENSG00000237613 0 0
ENSG00000268020 0 0
ENSG00000240361 0 0
ENSG00000186092 0 0
ENSG00000238009 17 2
ENSG00000233750 35 0
ENSG00000237683 8489 34
```

I checked the sample information file and there is no information about the conditions. I want to normalize the raw counts. For that, I want to use DESeq's `getVarianceStabilizedData()` function. However this function takes as an input a `CountDataSet` object. So when I try to make a `CountDataSet` object using this:

```
> cds <- newCountDataSet(countData = as.matrix(expr.df))
Error in is(conditions, "matrix") :
argument "conditions" is missing, with no default
```

It spits out an error asking me to specify the conditions. However, there are no conditions in this dataset. How can I normalize these values?

[1]: http://www.gtexportal.org/home/datasets

Answer: To save everyone's time, I tried this:

# specify just one condition
cds <- newCountDataSet(countData=as.matrix(expr.df),condition=c(rep(1,2921)))
cds <- estimateSizeFactors(cds)
# method ="blind" assumes replicates of a single condition
cds <- estimateDispersions(cds, method="blind", fitType="local")
vst <- getVarianceStabilizedData(cds) # this worked fine

Let me know if my approach seems reasonable or incorrect. I did not want to delete the question just in case other people face the same problem.
Question: **Tool used:** [seqkit][1]

Dummy fasta file (fasta.fa):

>test1
GCATCGATCAGCTACGATCATCACTA
GNNNNNNTACATCAGCACTACATCACTNNNNN
>test2
GTACGCTACGANNNGCTACGACTACGATATATATATATATATATATATATATATATATATATAT
GCTACGATCACNTACATCGACTA
>test3
GTGTGCTACATCATCACTACGTACTACAT
>test4
AA

**Command:**

./seqkit stat fasta.fa

**Output:**

file format type num_seqs sum_len min_len avg_len max_len
fasta.fa FASTA DNA 4 176 0 44 87

**Problem:**
**min_len** =0 (however, minimum length should be 2; sequence id "test4")

**Validation using seqkit:**

**Command:**

./seqkit fx2tab -l fasta.fa

**Output:**

test1	GCATCGATCAGCTACGATCATCACTAGNNNNNNTACATCAGCACTACATCACTNNNNN	58
test2	GTACGCTACGANNNGCTACGACTACGATATATATATATATATATATATATATATATATATATATGCTACGATCACNTACATCGACTA	87
test3	GTGTGCTACATCATCACTACGTACTACAT	29
test4	AA	2

**Notice: length of sequence test4 is "2"**

Is it a bug or I misunderstood something?

**PS: I am loving this tool (all thanks to [Wei Shen][2])and trying to exploit the utilities to make a new tool!**

[1]: https://github.com/shenwei356/seqkit
[2]: https://www.biostars.org/u/4664/

Answer: For the data you posted I get the correct result using Version: 0.3.4.1

file format type num_seqs sum_len min_len avg_len max_len
demo FASTA DNA 4 176 2 44 87
Question: I am using WGCNA for the first time to identify gene co-expression modules across a time course with 40+ samples (RNA-Seq). I've removed lowly expressed genes and focused the analysis on the 15,000 most dynamically expressed genes. The free-scale topology and other indicators (sample clustering, etc) look good and similar to other example datasets. The network is signed.

When I plot the "representative" merged eigengene module expression I get eigengene profiles for each module that agree well with my biological expectations. Nevertheless, when I look closer to investigate specific genes of interest in any particular module their expression pattern is completely discordant and sometimes entirely opposite to that of the "representative" eigengene.

Is this expected and if it is to what level? What could be causing this discordance? What parameters can I modify to make the coexpression clusters "tighter"? Is that necessary?

I suspect one of the factors that might be resulting is this pattern is eigenmodule merging. I used a merging distance threshold of 0.2 to merge eigenmodules but it might be a bit too strict although I'm not sure if there is a better way to choose this threshold. I am in fact expecting large numbers of modules as the conditions I'm comparing are fairly biologically different.

Find the module dendrogram below, as well as the merging cut-off (in red).

Any insights/pointers/suggestions would be greatly appreciated.

Thanks!

Answer: Some degree of discordance with a module's eigengene is expected. Have a look at the module membership values for these genes, it may well be that they're only borderline members. Regarding tightening up membership in a module, one simple method would be to cut the dendrogram at a lower value (exactly as you suspected). This will, of course, increase your overall module number, but since you're expecting that anyway it's unlikely to be an issue.

BTW, you can also look at some gene correlation heatmaps, which might help you determine a better cutting threshold.