cannin/biostars_qa huggingface dataset:
Question: Hello again Biostars community!

I am working on a project for my internship regarding the efficacy of different clustering techniques on RNA-seq data. I have used many different common clustering techniques in R and compared their results and the clusters they produce. Now, I have gotten to a point in the project where I want to analyze and explicitly show how my clustering is performing compared to "random" clustering. (Where I redistribute the "cluster" label across different groups of genes.)

To this point, I have identified k-means clustering as the most helpful for my dataset and identified some positive control groups that I know are biologically meaningful and should cluster together. For these positive controls, I have calculated a p-value and this is where my "random" groups come in!

For each of my clusters, I want to redistribute my genes as randomly as possible (within the limits of reproducibility) and examine how my p-values change for the positive control groups. Ideally I would like to see the p-values for my clusters showing much higher significance than the random clusters, but who knows? (I mean... I'm pretty sure I know, but that's why I'm doing this!)

I would like each new random cluster to contain the same number of genes as the clusters that I have already generated. (e.g. Cluster_1 had 60 genes, so Random_Cluster_1 should also have 60 genes etc.). Of course, I also need to make sure any one gene is only assigned to a single cluster.

Would anyone be able to recommend a robust method of "randomly" reassigning my genes to new clusters so that I can check my clustering performance?

Any feedback would be appreciated as I don't have much experience in this area (pretty new to bioinformatics!)

Thank you!

Answer: So you need to divide the genes (N=20000) into x clusters of different sizes. Each gene needs to be in one cluster.

Here some R solution (though not very optimal I guess - maybe there is already a function implemented in package that does that)

Example with N=20 genes and 4 clusters of size 2, 3, 5 and 10 genes.

assign_genes_randomly <- function(genes,cluster_sizes){
random_assignment <- sample(unlist(sapply(cluster_sizes,function(x){rep(x,x)})))
random_clusters <- sapply(cluster_sizes,function(x){genes[random_assignment==x]})
return(random_clusters)
}

genes <- 1:20
cluster_sizes <- c(2,3,5,10)

set.seed(123)
random_clusters <- assign_genes_randomly(genes,cluster_sizes)

Results :

random_clusters
[[1]]
[1] 6 20

[[2]]
[1] 12 18 19

[[3]]
[1] 1 3 9 11 14

[[4]]
[1] 2 4 5 7 8 10 13 15 16 17
Question: Hi,

I have a table like this;

6 6:29002062:rs7755402 0 29002062 G A
6 6:29004091:rs9468471 0 29004091 A G
6 6:29006250:rs9468473 0 29006250 A G
6 6:29006493:rs9461499 0 29006493 C A
6 6:29006844:rs7743837 0 29006844 G A

I want to remove everything before "rs" in the second column. I know I can use egrep,

egrep -o "(rs\S+)" file | cut -d " " -f 2 > newfile

However, then I`m left with only the rs string;

rs7755402
rs9468471
rs9468473
rs9461499
rs7743837
rs6919044
rs41424052
rs6924824
rs6456886
rs6456887

But I actually want the other columns too.
Any help is greatly appreciated!

Answer: Shell is great, but you will like [csvtk](https://github.com/shenwei356/csvtk), a cross-platform, efficient, practical and pretty CSV/TSV toolkit.

Using `csvtk replace` to **edit specific column(s)**, [download](http://bioinf.shenwei.me/csvtk/download/), [usage](http://bioinf.shenwei.me/csvtk/usage/#replace).

$ csvtk replace -H -t -f 2 -p '.+:' -r '' file
6 rs7755402 0 29002062 G A
6 rs9468471 0 29004091 A G
6 rs9468473 0 29006250 A G
6 rs9461499 0 29006493 C A
6 rs7743837 0 29006844 G A

long-option version:

$ csvtk replace --no-header-row --tabs --fields 2 --pattern '.+:' --replacement '' file
Question: Hi,
Is anyone aware of pubilcly available T Cell Receptor repertoire data sets from immuno-sequencing, particularly from tumor or PBMC in cancer patients or mouse. I am looking for a summarized form repertoire profile dataset that includes clonotype (AA and NT), number and frequency of clones, and V, D, J gene name for each clonotype.
Thank you,
- Pankaj

Answer: This [paper][1] may help.

[1]: http://www.ncbi.nlm.nih.gov/pubmed/26017500
Question: I have tried to break a problem down to two files. I want to take the most common 2nd column value for each sequence (1st column) then have the range of column 3 for those values. Also count the number of markers so the number of Fc06 for seq000001. Any idea how to do this? I was hoping for something easy on linux command line using awk but I just can't see it. The data is cM and chromosome file which I which to summarise.

file1
seq000001
seq000002

file2

seq000001	Fc06	35.948
seq000001	Fc06	36.793
seq000001	Fc06	37.009
seq000001	Fc06	37.009
seq000001	Fc06	37.009
seq000001	Fc06	37.009
seq000001	Fc06	37.009
seq000001	Fc06	37.009
seq000001	Fc06	37.009
seq000001	Fc06	37.009
seq000001	Fc06	37.009
seq000001	Fc06	37.009
seq000001	Fc06	37.009
seq000001	Fc08	37.009
seq000001	Fc06	37.368
seq000001	Fc06	37.368
seq000001	Fc06	37.368
seq000001	Fc06	37.368
seq000001	Fc06	37.368
seq000001	Fc06	37.58
seq000001	Fc06	37.916
seq000001	Fc06	37.916
seq000001	Fc06	37.916
seq000001	Fc06	37.916
seq000001	Fc06	37.916
seq000001	Fc08	37.916
seq000001	Fc06	37.916
seq000001	Fc06	37.916
seq000001	Fc06	37.916
seq000001	Fc06	37.916
seq000001	Fc06	37.916
seq000001	Fc06	38.1
seq000001	Fc06	38.664
seq000001	Fc06	38.923
seq000001	Fc06	38.978
seq000001	Fc06	39.324
seq000001	Fc06	40.341
seq000001	Fc06	40.341
seq000001	Fc06	40.341
seq000001	Fc06	40.543
seq000001	Fc06	40.543
seq000001	Fc06	40.543
seq000001	Fc06	40.992
seq000001	Fc06	41.189
seq000001	Fc06	41.189
seq000001	Fc06	41.347
seq000001	Fc06	41.39
seq000001	Fc06	41.39
seq000001	Fc06	41.39
seq000001	Fc06	41.399
seq000001	Fc06	41.399
seq000001	Fc06	41.399
seq000001	Fc06	41.657
seq000001	Fc06	41.71
seq000001	Fc06	41.785
seq000001	Fc06	41.923
seq000001	Fc06	42.237
seq000001	Fc06	42.634
seq000001	Fc06	42.963
seq000001	Fc06	43.285
seq000001	Fc06	43.478
seq000001	Fc06	43.478
seq000001	Fc06	43.478
seq000001	Fc06	43.478
seq000001	Fc06	43.744
seq000001	Fc06	43.744
seq000001	Fc06	45.234
seq000001	Fc06	45.234
seq000001	Fc06	46.272
seq000001	Fc06	46.272
seq000001	Fc06	46.272
seq000001	Fc06	46.272
seq000001	Fc06	46.272
seq000001	Fc06	51.086
seq000002	Fc06	63.754
seq000002	Fc06	63.754
seq000002	Fc09	13.078
seq000002	Fc09	13.078
seq000002	Fc09	13.078
seq000002	Fc09	16.342
seq000002	Fc09	16.342
seq000002	Fc09	16.342
seq000002	Fc09	16.342
seq000002	Fc09	17.004
seq000002	Fc09	17.004
seq000002	Fc09	17.004
seq000002	Fc09	17.004
seq000002	Fc09	17.004
seq000002	Fc09	17.004
seq000002	Fc09	17.004
seq000002	Fc09	17.004
seq000002	Fc09	17.004
seq000002	Fc09	17.004
seq000002	Fc09	17.806
seq000002	Fc09	18.544
seq000002	Fc09	19.59
seq000002	Fc09	19.59
seq000002	Fc09	19.59
seq000002	Fc09	19.59
seq000002	Fc09	19.59
seq000002	Fc09	19.59
seq000002	Fc09	19.59
seq000002	Fc09	19.59
seq000002	Fc09	19.59
seq000002	Fc09	19.59
seq000002	Fc11	19.59
seq000002	Fc09	19.59
seq000002	Fc09	39.857
seq000002	Fc09	20.558
seq000002	Fc09	20.892
seq000002	Fc09	21.323

Output

seq000001	Fc06	35.948-63.754 72
seq000002	Fc09	13.078-39.857 34

Answer: $ datamash -s -g 1,2 min 3 max 3
Question: I'm trying to align FASTA files with a custom genome reference (GTF file), which contains a few sequences for genes that were inserted. So far I've just run cellranger count on these files and everything exited without errors. However, in the downstream analysis the result is not as expected, and I am wondering whether this could have to do with the alignment.
Basically, I am detecting the inserted sequences in many cells that should not contain it. One possibility is of course experimental error or contamination, but I was thinking it could also have to do with the alignment misaligning reads for those inserted sequences and detecting counts where there should not be counts.

Hence I was wondering **if there are parameters to tune the "sensitivity" of the genome alignment** in cellranger. Basically, I would want to be strict and discard reads which have even a small number of mismatches. Or if you have other suggestions to filter out potential misaligned reads, I'd be happy to know as well.

Cellranger count itself does not seem to have any relevant options, but it seems that on the background it runs STAR. Hence if it is possible to pass STAR arguments to cellranger, that could also be a possible solution (although I haven't found the relevant STAR parameters yet). Or maybe I'm just misunderstanding genome alignment altogether and there is not much to tune with the current algorithms?

Answer: There is no way to do this via cellranger without [altering the source code](https://kb.10xgenomics.com/hc/en-us/articles/360003877352-How-can-I-modify-the-STAR-alignment-parameters-in-Cell-Ranger-). You could always trying other gene quantification/alignment options though, like [alevin](https://salmon.readthedocs.io/en/latest/alevin.html) or [STARsolo](https://github.com/alexdobin/STAR/blob/master/docs/STARsolo.md) which would allow you full control over the parameters used. STARsolo even has directions for how to largely replicate cellranger results, so you could do that and tweak only your parameters of interest.
Question: I have a question about how to change order of clusters on heatmap.
I tried to make heatmap using Complex Heatmap package and I made K-mean clusters. my code is below:

HM <- Heatmap(matrix2, km = 4, show_row_names = FALSE, show_column_names = F, top_annotation = ha_column, cluster_columns = FALSE, show_heatmap_legend = T)
draw(HM, heatmap_legend_side = "right", show_annotation_legend = F)

I got this heatmap as below.
https://ibb.co/jqGb77

In this example, I want to move cluster 2 to top, cluster 1 in 2nd top, cluster 3 on 3rd and cluster 4 on the bottom. How can I re-order the position of clusters ?

Answer: I'm a big fan of *ComplexHeatmap* and again kudos to the developer for making it such a flexible package. However, I'm not immediately aware of a function that does (directly) what you want.

I would prefer to perform the k-means outside of the *ComplexHeatmap* package by just using `kmeans()`, and then manipulating the gene-to-cluster assignment vector in order to dictate the order of the clusters.

Create some random data (20 genes x 20 samples):

test <- matrix(rexp(400, rate=.3), ncol=20)
rownames(test) <- paste(rep("Gene", nrow(test)), 1:nrow(test), sep="")
colnames(test) <- paste(rep("Sample", ncol(test)), 1:ncol(test), sep="")

We then perform k-means with 4 centers and take a look at the gene-to-cluster assignment:

kclus <- kmeans(test, 4)
kclus$cluster
Gene1 Gene2 Gene3 Gene4 Gene5 Gene6 Gene7 Gene8 Gene9 Gene10 Gene11
3 3 2 1 2 1 3 1 2 1 4
Gene12 Gene13 Gene14 Gene15 Gene16 Gene17 Gene18 Gene19 Gene20
4 4 1 2 4 2 2 4 2

We then use the *split* parameter of `Heatmap()` in order to split the heatmap based on the k-means result:

split <- paste0("Cluster\n", kclus$cluster)
default.hmap <- Heatmap(test, split=split, cluster_row_slices = FALSE)

As we want to fix the order of the clusters, we have to re-order the gene-to-cluster assignment as a factor:

split <- factor(paste0("Cluster\n", kclus$cluster), levels=c("Cluster\n2","Cluster\n1","Cluster\n3","Cluster\n4"))
reorder.hmap <- Heatmap(test, split=split, cluster_row_slices = FALSE)

Then draw both the default-ordered and then re-ordered:

pushViewport(viewport(layout=grid.layout(nr=1, nc=2)))

pushViewport(viewport(layout.pos.row=1, layout.pos.col=1))
draw(default.hmap, newpage=FALSE)
upViewport()

pushViewport(viewport(layout.pos.row=1, layout.pos.col=2))
draw(reorder.hmap, newpage=FALSE)
upViewport()

upViewport()
Question: Hello,
I am trying to extract reads at 1 base pair from many bam files. Each bam file has a bai file with the same name+.bai in the same folder.

When I run

samtools view myfile.bam

I get the output, but when I add coordinates such as

samtools view myfile.bam "chr1:1-100000"

I get the following error:

[bam_index_load] fail to load BAM index.
[main_samview] random alignment retrieval only works for indexed BAM files.

My understanding is that the bai file is the index it needs. Am I missing something?

Thanks in advance.

Answer: Some things to consider:

1. Your index may old or generated with a program other than samtools. You might want to regenerate it.
2. Make sure the BAM itself is sorted by coordinate. Also make sure that the names are correct between the two files, i.e., you didn't change the name of the BAM after indexing or something.
3. The index itself also needs to be accessible and readable by samtools.
Question: I am conducting differential expression analysis in limma for a time course microarray experiment. I have 3 time points (1,2,3) and two conditions (control,disease). My contrast is formulated as:

```
pheno<-factor(c("disease.1","disease.1","disease.2","disease.2","disease.3","disease.3","control.1","control.1","control.2","control.2","control.3","control.3")
​model<-model.matrix(~0+pheno)
```

```
contrast<-makeContrasts(disease.1-control.1, disease.2-control.2, disease.3-control.3,levels=model)
fit<-lmFit(expression,model)
fit2<-contrasts.fit(fit,contrast)
fit2<-eBayes(fit2)

topTableF(fit2)
```

How can I apply SVA with this experimental setup? Or is there an alternative? The challenge I am facing is SVA seems designed for two class experiments, because the null model requires leaving out the variable of interest, however in this case there isn't one variable but many.

Thanks!

Answer: You apply the sva function to a normalised matrix of expression values. I'm assuming that you're just looking for potential surrogate variables and not trying to correct for a known batch effect?

Edit: Also, I'd recommend you check out the time course section of [the limma user guide][1]

[1]: http://www.bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf
Question: Hello friends, please help me

i entered these codes, but whenever i use the output as input to geworkbench, it says that" could not parse line #1, line should have 12 columns but has 13 columns, which need manually add a tab at the beggining of the header to make it a valid RMA format"

library(affy)

> data <- ReadAffy()
> exp <- rma(data)

Background correcting
Normalizing
Calculating Expression
> control <- grep("AFFX",rownames(exp))
> exp <- exp[-control,]
> e <- exprs(exp)
> e <- apply(e,1,function(x) x-median(x))
> e.scaled <- t(scale(t(e),center=FALSE))
> write.exprs(exp, file="output.txt", sep="\t")

where i did wrong please??

Answer: Change:

write.exprs(exp, file="output.txt", sep="\t")

to:

write.exprs(exp, file="output.txt", sep="\t", row.names=FALSE)
Question: I have multiple VCFs for single cell RNA seq data and I want to get peptide sequences from these files. I have searched and found out that we can annotate the VCF with Ensembl's VEP to get the amino acids/protein information. However, I am looking to get a fasta file as an output which can be used for downstream analysis. I found out that GATK's FastaAlternateReferenceMaker can be used to get a fasta file from VCF.

Can I use the output VCF from VEP as an input to the FastaAlternateReferenceMaker to get the required fasta? I am not sure if I should pass the entire VCF as an input or just the protein information. Please help me get a better understanding of this procedure.

Answer: I found out that this can be done in the following way:

1. Find the proteins in the mutations using Ensembl's VEP and sending the VCF files as an input. By doing so we will get several information as an output but we are interested in Ensembl’s VEP protein ID which starts with **ENSP.** (VEP output format provided here: [https://m.ensembl.org/info/docs/tools/vep/vep_formats.html#vcfout][1])
2. Once we have the protein IDs we can use the REST API provided by Ensembl and pass the ID to get the fasta:
[https://rest.ensembl.org/sequence/id/ENSP00000404426?content-type=text/x-fasta;type=protein][2]

[1]: https://m.ensembl.org/info/docs/tools/vep/vep_formats.html#vcfout
[2]: https://rest.ensembl.org/sequence/id/ENSP00000404426?content-type=text/x-fasta;type=protein
Question: I want to extract a list of sequences from NCBI. I am esearch command mentioned [here][1]. For one gene symbol, I could do it like this:

esearch -db nuccore -q 'SS1G_01676[gene]' | efilter -source refseq -molecule genomic | efetch -format gene_fasta | awk -v RS='(^|\n)>' '/SS1G_01676/{print RT $0}'

I want to use a bash loop to extract a list of sequences and below is what I have tried, but wouldn't yield any results. What am I missing here?

declare -a arr=("SS1G_03709" "SS1G_07286" "SS1G_04907")
for i in "${arr[@]}"
do
myquery="'${i}[gene]'"
echo "myid :" ${i}
echo "my query :" ${myquery}
esearch -db nuccore -q ${myquery} | efilter -source refseq -molecule genomic | efetch -format gene_fasta | awk -v RS='(^|\n)>' '/${i}/{print RT $0}' >>text.fasta
done

[1]: https://www.biostars.org/p/345772/#345918

Answer: You are right, it stops after the first gene. I am not sure why. (EDIT: see explanation below) This works though:

for line in `cat temp.txt`; do esearch -db nuccore -q $line| efilter -source refseq -molecule genomic | efetch -format gene_fasta | awk -v r="$line" 'BEGIN {RS="(^|\n)>"} $0 ~ r {print ">" $0}' ; done

Apparently, for while loops involving `esearch` one should add `< /dev/null` to prevent esearch from "draining" the remaining lines from stdin. See documentation here: https://www.ncbi.nlm.nih.gov/books/NBK179288/#chapter6.While_Loop
Question: Hi,

I have a list of kmers, between 8-12 nt in length, and I would like to align these to a larger sequence returning all ungapped matches with at most 2 mismatches. I would like search to be exhaustive i.e. I do not want to miss anything. I wrote a python script to compute hamming distance for all substrings of my reference to the query, but it is too slow for many (1000s) queries on a reference of ~100,000nt.

What program would you recommend that does this and runs rather quickly. I have looked into Bowtie2, but I am unsure if it was designed to work with such short query sequences.

Thanks for the feedback.

Answer: Ended up creating a custom python script that:

1. Broke up the larger sequence into kmers of required sizes into a set.
2. For each query kmer, compute possible all 2-mismatch kmers into a set.
3. interesect set from 1 with set from 2.

Works extremely quickly as my kmers are quite small (8-12 nt) and the search target is also relatively small (tens of kb).
Question: Hi all,

There is a great preprint for the "complete" human genome T2T-CHM13 on bioRxiv. I'm not able to find a gff/gtf/bed annotation file for this and I was hoping someone might share a link if one exists.

[https://www.biorxiv.org/content/10.1101/2021.05.26.445798v1.full][1]

[https://github.com/marbl/CHM13][2]

Thanks!

[1]: https://www.biorxiv.org/content/10.1101/2021.05.26.445798v1.full
[2]: https://github.com/marbl/CHM13

Answer: Hey, it looks like the correct link is now on Github.

See
https://github.com/marbl/CHM13#downloads.
Question: Hi Everyone,

I am a newbie in processing scRNA-Seq data and I am trying to understand a Seurat Object.
Here's a single cell dataset deposited on NCBI GEO: GSE158130

```
str(obj)

Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
..@ i : int [1:2110714] 25 26 28 30 32 112 113 114 124 127 ...
..@ p : int [1:1177] 0 2014 4094 6011 8780 11407 14821 17657 19653 20675 ...
..@ Dim : int [1:2] 27044 1176
..@ Dimnames:List of 2
.. ..$ : chr [1:27044] "1/2-SBSRNA4" "5S_rRNA" "5_8S_rRNA" "7M1-2;OR2F1" ...
.. ..$ : chr [1:1176] "WSJ0005001" "WSJ0005002" "WSJ0005003" "WSJ0005004" ...
..@ x : num [1:2110714] 1 1 2 1 1 1 1 1 1 1 ...
..@ factors : list()
```

What I understand from this object:
- @Dim provides the dimensionality of this dataset is 27,044 (genes) x 1176 (cells)
- @Dimnames slot provides names for genes and cells.

I have the following questions:
1. What data does slot @i, @p and @x show?
2. Could we determine if QC, normalization and dimensionality reduction methods have been applied to the data by looking at the object?
3. This dataset also provides 2 additional files - `GSE158130_cellID_barcode_map.txt` and `GSE158130_SK-N-SH_counts.txt`. I am not sure how to use these files and what should be my next steps.

Any help is appreciated here.

Thanks.

Answer: This isn't a Seurat object, but rather a sparse matrix from the `Matrix` library. `@x` are all of the non-zero counts. `@i` are the index positions column-wise of non-zero values. `@p` is a cumulative sum of the number of non-zero values in each column.

`GSE158130_SK-N-SH_counts.txt` is the same matrix but in csv format. `GSE158130_cellID_barcode_map.txt` has some metadata about each cell, such as the assigned cell-id, corresponding cell barcode, patient ID, etc.

From this data alone you can't tell how the data was processed. You'll need to look at the methods from the paper.
Question: Hi everyone.

I am having troubles undestanding if the PFAM-A hmm's provided by the pfam database are the ones obtained by the SEED or from the FULL alignment of the family.

In the step-by-step description in pfam's website (say: http://pfam.xfam.org/family/PF00004#tabview=tab6 ) they provide the command used to obtain the hmm form the seed, and then the command to run the 'search' performed on all the pfamseq to obtain the "full" msa. Then they provide a "raw hmm" to download, but it is not clear to me if this hmm is the one they obtained from the seed with the given command OR another one, obtained from the full alignment; neither I unsterstand if it is supposed to be the same hmm i can find in the ftp website in the file PFAM-A.hmm

In the PFAM papers i cannot find any paragraph citing the "full-msa" hmm, but I've head about it from many people. Maybe they are wrong Or i am not good in reading pfam's papers and documentation, but it's days i am struggling with this and I still cannot find some clear definition of all the databases and the hmms used to build them.

Answer: I was once confused by this. It seems definitions are not as clear in bioinformatics as they are in mathematics and misunderstanding is common. If you click on the 'download' link once you follow your link above you will see the NSEQ field which has the value 207. This seed alignment defines the family.

Sometimes the HMM model is called the alignment, sometimes all the sequences lined up together is called the alignment. The full alignment is what results when executing an HMM search based on the seed alignment. The full alignment does not define a family since it changes as new sequence data is ingested.
Question: Hi folks. I need to run a de novo short-read genome assembler (on a paired-end/mate-pair library) that prefers outputting shorter but error-free contigs rather than longer contigs/scaffolds which may be mis-assembled. What assembler or what specific setting in an assembler of choice do you recommend to yield such contigs (as error-free as possible and no contig overlappings)?

Answer: I think error free contigs depends on the quality of your data too and the contamination if any. It also depends on the repetitiveness of genome, level of polymorphism (in order to know the correctness of contigs) and heterozygosity of the individual. SOAP contigs are short as they start from K+1 of your kmer. By increasing the `min_abundance` parameter in denovo assemblers, you can get more accurate contigs. Minia is definitely one of the ones to try out.

If you have lesser number of error-free reads, go for overlap assembler such as CAP3. This wouldn't work for a large number of reads due to memory constraints.