text:
if, else statment gives back else value
I have a list of protein sequence ids (names). I want to sort the lines of a table, checking if the id is equal to the first column of the table. I split the table on lines and search for the id via regex:
code:
def map_line_to_ids(id):
    id= id.strip('\n')
    for s in lines:
        if re.search(id,s):
            return s #<------------ why always return else not if(s) when takes empty 'def'
        else:
            return id + '\n'

names = ['XXXXXXX45.1', 'YP_343434.1', 'WP_XXXXXXXX2.1', 'WP_XXXXXXX90.1',
         'WP_0XXXXXX62.1', 'XP_222222222.2', 'WP_192XXXXXX.1',
         'WP_XXXXXXX29.1', 'WP_0XXXXXXXX.1']

sorted_lines = list(map(map_line_to_ids, names))

names (order):
XXXXXXX45.1
YP_343434.1
WP_XXXXXXXX2.1
WP_XXXXXXX90.1
WP_0XXXXXX62.1
XP_222222222.2
WP_192XXXXXX.1
WP_XXXXXXX29.1
WP_0XXXXXXXX.1

Input table:
WP_XXXXXXXX2.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
WP_192XXXXXX.1  PF13489 Methyltransferase domain    4
XXXXXXX45.1 PF02310 B12 binding domain  12; PF04055 Radical SAM superfamily 5
WP_XXXXXXX29.1  PF13489 Methyltransferase domain    4
WP_0XXXXXXXX.1  PF13489 Methyltransferase domain    4
WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5; PF06325  Ribosomal protein L11 methyltransferase (PrmA)  7
WP_0XXXXXX62.1  PF13489 Methyltransferase domain    4

Expected output:
XXXXXXX45.1 PF02310 B12 binding domain  12; PF04055 Radical SAM superfamily 5
YP_343434.1
WP_XXXXXXXX2.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5; PF06325  Ribosomal protein L11 methyltransferase (PrmA)  7
XP_222222222.2
WP_0XXXXXX62.1  PF13489 Methyltransferase domain    4
WP_192XXXXXX.1  PF13489 Methyltransferase domain    4
WP_XXXXXXX29.1  PF13489 Methyltransferase domain    4
WP_0XXXXXXXX.1  PF13489 Methyltransferase domain    4

Hi @MTG, this is a pure programming question, so it would have been appropriate to ask on StackOverflow. But I think we're happy to handle it here on Bioinformatics as well.

Following on your last post, your initial file test.txt looks like that:
WP_XXXXXXXX2.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
WP_192XXXXXX.1  PF13489 Methyltransferase domain    4
XXXXXXX45.1 PF02310 B12 binding domain  12
XXXXXXX45.1 PF04055 Radical SAM superfamily 5
WP_XXXXXXX29.1  PF13489 Methyltransferase domain    4
WP_0XXXXXXXX.1  PF13489 Methyltransferase domain    4
WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  7
WP_0XXXXXX62.1  PF13489 Methyltransferase domain    4

And you have another file order.txt that contains your ordered IDs:
XXXXXXX45.1
YP_343434.1
WP_XXXXXXXX2.1
WP_XXXXXXX90.1
WP_0XXXXXX62.1
XP_222222222.2
WP_192XXXXXX.1
WP_XXXXXXX29.1
WP_0XXXXXXXX.1

Since you want to merge by the first field in test.txt and sort them according to order.txt, you have to read first both of the files:
# read 'test.txt'
with open("test.txt",'r') as rf:
    l = rf.readlines()

# read sorted IDs 'order.txt'
with open("order.txt",'r') as rf:
    l_order = rf.readlines()

# remove newlines from each element in list
l_order = [i.strip() for i in l_order]

My previous answers describe how we obtain the merged fields in the file, so I will not go over it again. Since we need to sort the dictionary according to the list, we simply have to do:
# for each element in ordered list
for i in l_order:
    if i in dict_id.keys():                        # if element is in my dictionary's keys
        print(i + "\t" + "; ".join(dict_id[i]))    # print the key and the value
    else:
        print(i)                                   # else prints the element

This should print:
XXXXXXX45.1     PF02310 B12 binding domain      12; PF04055     Radical SAM superfamily 5
YP_343434.1
WP_XXXXXXXX2.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5; PF06325      Ribosomal protein L11 methyltransferase (PrmA)  7
WP_0XXXXXX62.1  PF13489 Methyltransferase domain        4
XP_222222222.2
WP_192XXXXXX.1  PF13489 Methyltransferase domain        4
WP_XXXXXXX29.1  PF13489 Methyltransferase domain        4
WP_0XXXXXXXX.1  PF13489 Methyltransferase domain        4

Full code:
with open("test.txt",'r') as rf:
    l = rf.readlines()

with open("order.txt",'r') as rf:
    l_order = rf.readlines()

l_order = [i.strip() for i in l_order]

dict_id = {}
for i in l:
    key_id = i.split("\t")[0]
    if key_id not in dict_id.keys():
        value = "\t".join(i.split("\t")[1:]).rstrip()
        dict_id[key_id] = [value]
    else:
        dict_id[key_id].append("\t".join(i.split("\t")[1:]).rstrip())

for i in l_order:
    if i in dict_id.keys():
        print(i + "\t" + "; ".join(dict_id[i]))
    else:
        print(i)

PS: please for further pure programming questions, do post your questions in StackOverflow instead of Bioinformatics.


text:
Integrating scRNA-seq data using raw data
I am trying to reporduce the the following article  Chen, M. B., Jiang, X., Quake, S. R., & Südhof, T. C. (2020). Persistent transcriptional programmes are associated with remote memory. Nature, 587(7834), 437-442..  In the methods section Bioinformatics and data analysis - Mapping to the genome  they say "Next, the ‘canonical correlation analysis’ function from the Seurat package27 was used to align raw data from multiple experiments. "  But they had previously aligned the data and generated gene counts.  So by this do they mean anchor?  I'm confused why would some one run CCA on the raw unaligned data.  They only provided gene count data not raw data, but they use the CCA to restrict the gene's used from ~23000 down to just ~2500 of the more positively and negatively correlated (to the resulting CCA generated principle components) without providing the exact details of which genes they included or excluded before running their analysis.
Can I just do the CCA on the count data and hope for the same results?  Have I misunderstood their filtering and alignment?
EDIT  For clarity for future readers, my question boiled down to the the term alignment. I was confused because I thought alignment specifically referred to aligning sequenced fragments to a reference genome, here alignment basically refers to correcting for batch effects in scRNA-seq.

I believe when you say alignment, you mean aligning reads to a genome (sometimes to transcriptome) and count these to get count matrices. In the aforementioned paper, however, what is meant is "bringing different data sets to a level where they can be compared/integrated/...". Basically scRNA-seq data are heavily prone to batch effetcs and if you try to integrate different data sets (let's say data generated with 10x 3' and 10x 5' chemistries), there is a good chance is most of the signal would come from the technical differences rather than biology, for example resulting in clustering by scRNA-seq technology and not cell type/stat, hence there are many batch-effect-removal/data-integration/data-alignment tools for scRNA-seq.
Seurat has a nice CCA tutorial, you would not need the raw data, count data would be fine but you will need to be able to specify what data are coming from what source, the authors must have provided this information.
When it comes to what an "anchor" is, you would need to read about the CCA technique but in a rather naïve way, you can think of these as commonalities between data sets or some kind of a reference and are used to "align" these different data sets.
Moreover, reducing the number of genes from tens of thousands to thousands is just fine and is not really specific to CCA: not all genes are expected to be expressed at all times in all cell types and hence are not that informative and also scRNA-seq is not able to pick up lowly-expressed genes, causing these to show up with 0 counts in the count matrix, so again not that informative/interesting for the downstream analysis steps.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
How to get the number of complete phage genomes available on ncbi?
I am looking to establish the total number of complete phage genomes available on NCBI.
I am not looking for any specific type, but want to understand the total diversity available.
How can I perform the search?

Could you supply more details? What groups of phages? Any in particular? Every phage? The answer depends on what you are looking for.

The idea is to get the total number of complete phage genomes available on ncbi. So yes, I am not looking for any particular kind, just want to know how many are there.

That's easy. Please update your question to clarify this point, because its helps the archives.

Done, thank you for guiding me to make the question better suited for the forum.

Go to the NCBI nuccore database, at https://www.ncbi.nlm.nih.gov/nuccore
Copy and paste the search term:
(bacterial virus[Organism]) AND (3000:110000[sequence length]) NOT bacteria[Organism]

You'll get 74583 results.

Click on the link send to and click the checkbox called file.
It will then download.

Check the downloaded file is in fasta format with ">" to start the ID.

If fasta format is absent, return to the webpage and where it states summary switch that to fasta text

I can't remember if point 4 is needed.
Notes

Sequence length, I am guessing that all phages are between 3000 and 105000 in length, this is to screen out partial sequences and to minimise  bacteria creeping in to the search. I don't quite understand how a phage can be 105 000bp in length, but you can twiddle these parameters to refine the search.
@acvill informs us that giant phages greater than >200 kb have been annoated. Thus a separate search is likely required targeting this genome size.

Anecdotally, (3000:110000[sequence length]) should capture the vast majority of phage genomes, though a growing number of huge phages (>200 kb) are being annotated.

Thank you, that is good expert opinion. I will update the post. I am mystified how a virus can be this size in comparison to their bacterial host genome, but phages are not my thing.

I did not know you could download all the search results at once as fasta.

How about trying this?

Go to https://www.ncbi.nlm.nih.gov/nucleotide. Then type the following command in the search box. It resulted in 23971 items.

(bacterial virus[Organism]) AND (complete genome[ti]) NOT shotgun[ti] NOT plasmid[ti] NOT bacteria[Organism]

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Filtering a vcf with a text file of SNP rsIDs
I have a vcf file containing an ID field with a snps and a list of human SNPs from dbSNP.I want to remove whatever in dbSNP from vcf file. The text file is like this:
rs978760828
rs1008829651
rs1052373574
rs1326880612
rs768019142
rs892501864
rs1010989343
rs1489251879
rs1022805358
rs1322538365

How to remove the common snps from dbSNP and my vcf?

I'd advise against using grep, sed, awk for processing vcfs unless you really know what you're doing - there's almost always an easier way using bcftools.
bcftools view -i'ID=@snplist.txt' in.vcf


text:
Separation of dihedral angles of the desired sequence in the protein
have a list of pdbid and I want to separate the dihedral angles of some sequences in protein .
pdbid      st          en
101m    4.A.GLU 15.A.ALA
102m    4.A.GLU 15.A.ALA
103m    4.A.GLU 15.A.ALA
104m    4.A.GLU 15.A.ALA
105m    4.A.GLU 15.A.ALA
106m    4.A.GLU 15.A.ALA
107m    4.A.GLU 15.A.ALA

st is start of sequences   en is end of sequences.

And after separating the dihedral angles, I want to write the
corresponding pdbid in front of it.
I wrote the following code:
  library(bio3d)

  # PDB structure object as obtained from function read.pdb.
 pdb=read.pdb(df$pdbid[1])

  #Calculate all torsion angles for a given protein PDB structure object
tor=torsion.pdb(pdb)

            phi        psi
1.A.VAL -144.91395  137.473333
2.A.LEU -95.02097   149.7834
3.A.SER -82.68263   166.676942
4.A.GLU -62.53142   -36.263406
5.A.GLY -59.64632   -42.035533
6.A.GLU -68.69144   -40.23538
7.A.TRP -63.52969   -36.769049
8.A.GLN -64.80557   -38.512821
9.A.LEU -64.82472   -42.143097
10.A.VAL    -63.77383   -50.465974
11.A.LEU    -69.74235   -25.698665
12.A.HIS    -70.86863   -41.810178
13.A.VAL    -74.99987   -33.340633
14.A.TRP    -63.20999   -32.690374
15.A.ALA    -61.13164   -30.030673
16.A.LYS    -75.94714   -43.271738
17.A.VAL    -54.76125   -43.181171
18.A.GLU    -64.33161   -12.649873
19.A.ALA    -70.38802   -22.502227

#separating phi angles from "4.A.GLU" to "15.A.ALA"
phi=tor$phi[df$st[1]:df$en[1]]

#separating psi angles
psi=tor$psi[df$st[1]:df$en[1]]
print(cbind(phi,psi))

#Put the codes in the loop
for(i in 2:nrow(df))
{
  pdb=read.pdb(df$pdbid[i])
  tor=torsion.pdb(pdb)
  psi=tor$psi[df$st[i]:df$en[i]]
  phi=tor$phi[df$st[i]:df$en[i]]
  print(cbind(phi,psi))
}

But I get the following error:
Error in df$st[1]:df$en[1] : NA/NaN argument
In addition: Warning messages:
1: NAs introduced by coercion
2: NAs introduced by coercion

str(df)
'data.frame':   9 obs. of  3 variables:
 $ pdbid: chr  "101m" "102m" "103m" "104m" ...
 $ st   : chr  "4.A.GLU" "4.A.GLU" "4.A.GLU" "4.A.GLU" ...
 $ en   : chr  "15.A.ALA" "15.A.ALA" "15.A.ALA" "15.A.ALA" ...
str(tor)
List of 11
 $ psi   : num [1:295] 158.1 137.5 149.8 166.7 -36.3 ...
 $ phi   : num [1:295] NA -144.9 -95 -82.7 -62.5 ...
 $ omega : num [1:295] -179 179 177 179 179 ...
 $ chi1  : num [1:295] 54.3 175.7 -81 62.3 175.1 ...
 $ chi2  : num [1:295] -171.5 NA 72.7 NA 64.4 ...
 $ chi3  : num [1:295] -72.1 NA NA NA 17 ...
 $ chi4  : num [1:295] NA NA NA NA NA NA NA NA NA NA ...
 $ chi5  : logi [1:295] NA NA NA NA NA NA ...
 $ alpha : num [1:295] NA -139.4 -123.2 -108.3 58.3 ...
 $ coords: num [1:3, 1:17, 1:295] 24.28 8.37 -9.85 24.4 9.86 ...
  ..- attr(*, "dimnames")=List of 3
  .. ..$ xyz: chr [1:3] "x" "y" "z"
  .. ..$ atm: chr [1:17] "N" "CA" "C" "O" ...
  .. ..$ res: chr [1:295] "  0.A.MET" "  1.A.VAL" "  2.A.LEU" "  3.A.SER" ...
 $ tbl   : num [1:295, 1:7] NA -144.9 -95 -82.7 -62.5 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:295] "  0.A.MET" "  1.A.VAL" "  2.A.LEU" "  3.A.SER" ...
  .. ..$ : chr [1:7] "phi" "psi" "chi1" "chi2" ...

Welcome to the site. This is a good report that made it easy to spot the error (the only missing piece is what did you try to solve this). There seems to be a problem on tor$phi[df$st[1]:df$en[1]] . Looking at the sample data they are characters, but you are using them as numbers, I'm not sure if the packages provides some kind of function to use them like that, but my guess is that it doesn't.

NAN : not a number  ? but NA ? is it null value ?

have you tried asking here : https://bitbucket.org/Grantlab/bio3d/issues


text:
Translating Simpson and Shannon diversity indices?
So I am trying to compute diversity indices for a dataset I have but the math doesn't work out and I don't really understand how to translate Shannon and/or Simpson diversity indices into something that is understandable. Maybe I got everything wrong, this is certainly a possibility. However, this is what I have up to now:
# load libraries
library("phyloseq")
data("GlobalPatterns")

# load data
globp = GlobalPatterns

# compute transformation according to link below
simpson_transformed <- 1/(1-estimate_richness(globp, measures="Simpson"))
shannon_transformed <- exp(estimate_richness(globp, measures="Shannon"))

# print side by side
print(cbind.data.frame(simpson_transformed, shannon_transformed))

The comparing_indices columns are computed according to what I found at this link (italic is mine):

They [i.e. the values] should be converted to effective number of species,
which are the true diversities. By referring to Table 1 we see that
species richness is already a true diversity, so the true diversity
according to this index is 5.000 species. Table 1 shows that the
Shannon entropy or Shannon-Wiener index is converted to effective
number of species or true diversity by taking the exponential:
exp(1.609) is 5.000, so the true diversity according to the Shannon
entropy is also 5.000 species. The Gini-Simpson index is converted to
a  true diversity by subtracting it from unity and inverting:
1/(1-0.8) = 5.000 species also. So in fact all these indices agree
that the diversity of this community is 5.000 species.

I also verified how phyloseq is computing the indices. I believe the transformation for Shannon should be the same mentioned above. I have doubts about the Simpson, since they refer to the Gini-Simpson index and not the Simpson itself. However, the documentation for vegan's diversity() function, which is the function that (I believe) is called by phyloseq's estimate_richness() (see docs), mention that the Simpson index is returned as 1-D which seems to be the Gini-Simpson index (see the 1-lambda). In other words, it looks like the transformation I linked above should be exactly the one that are used in the code I provided. But, I was expecting the same values in the two columns. But they are not there and I wonder why this happens. What did I do wrong?

You are trying to obtain Hill numbers, sometimes called true diversity, and you have different results because you have uneven species frequencies.
Vegan calculates Shannon index like:
$$ H' = -\sum_{i=1}^R p_i \ln p_i $$
And Gini-Simpson like:
$$ 1 - \lambda = 1 - \sum_{i=1}^R p_i^2 $$
So the Hill numbers should be:
$$ H1 = exp(H') $$
$$ H2 = 1 / 1 - \lambda $$
With Phyloseq:
H1 <- 1 / estimate_richness(globp, measures="Simpson")
H2 <- exp(estimate_richness(globp, measures="Shannon")

I do have a question, though...shouldn't H1 and H2 be equal to each other, since we are talking about true diversity? Why is that not the case?

H1 is equal to H2 for an even species frequency. The text that you cite uses the term "true diversity" or the "effective number of species", which means Hill numbers.

alright I think I am understanding what you meant. thank you!


text:
How to manually curate a genome assembly for sequence variation or error?
I have a PacBio HiFi assembly of 1.1 Gb from a heterozygous species. I have aligned this assembly against a reference genome which is around 0.9 Gb. I can see that there are quite a few INDELs, Duplications, Inverted regions in my Query sequence on the Y-axis (PacBio Hifi) compared to the reference genome.
I want to check if these structural variations are real or only resulting from sequencing errors or not. So, I want to manually curate the PacBio assembly. Because I can see at which contig the structural variations are.
Can you please suggest any tools for manual curation of a genome assembly to investigate whether these structural variations are sequencing errors or real sequence variation?
I have uploaded the dot-plot below. Thank you.

based on your other question, I would suggest mapping the hi-c reads to this assembly and manually polishing in juicebox

So, you mean take the HiFi fasta and map the Hi-C fasta on it? Or do I need raw reads from the Hi-C reads?

You would need to map the raw Hi-C reads. the FASTA is the outcome of a lot of procedures that remove information. It is the Hi-C reads have the actual proximity information about which regions are close to each other on the chromosome. You can use that to correct structural errors (if any). If those are truly errors then they will look weird on the Hi-C heatmap and you can correct them manually in juicebox: https://aidenlab.gitbook.io/juicebox/assembly-tools

Ah ok. Thanks for clarifying the confusion. It is clear now.


text:
Is there a way to get a list of all PDB IDs for proteins with crystal structure information?
I would like to download all the PDB files which have crystal structure information in them. But the query builder seems to require me to enter the PDB IDs manually.
Is there a way to get a list of those PDB IDs at RCSB PDB?
Thanks in advance.

Thank you for the links. Very useful indeed!

Comment migrated as answer as it inadvertently may have been an answer.

If the sole records in a PDB file that are of interest are EXPDTA,CRYST1 ORIGX1``,ORIGX2, ORIGX3,SCALE1,SCALE2,SCALE3, the in the FTP site the file ftp.wwpdb.org/pub/pdb/derived_data/index/crystal.idx  contains this information in a tabular format.
If an extra piece of metadata is sought, such as some of the refinement details in REMARK lines, then the PDBe API can help. In particular, the experiment entrypoint (example) is of interest.

running the query builder at RCS PDB against id with space group (so crystal info) I get:

more about RCSB PDB web services here: https://www.rcsb.org/docs/programmatic-access/web-services-overview
and the search API here: https://search.rcsb.org/index.html#search-api  think you should query for entries where symmetry.pdbx_full_space_group_name_H_M exists (https://search.rcsb.org/structure-search-attributes.html and ask to have returned entry.id (a list of in this case)


text:
tensorflow nn_model for DNA sequences: Matrix size-incompatible: In[0]: [2,1], In[1]: [784,300]
Hope anyone can help a beginner here. I'm building a proof-of concept tensorflow classifier for DNA sequences. However, the NN model does not let through train and test vectors saying the matrix size is incompatible. Where could be the mistake in my implementation?
import numpy as np
import pandas as pd

df=pd.DataFrame(columns=['sequences', 'label'])
a='ttttccagaattctcttagttt gtgatgtctttattgcttctattt'
b='ctcctgcttgctttttttcttg ggtttctgatattctttaaaggat'
c='tcctgcttgctttttttcttgg gtttctgatattctttaaaggatt'
df.sequences=[a,b,c]
df.label=[1,1,1]

from sklearn.preprocessing import LabelEncoder, OneHotEncoder
integer_encoder = LabelEncoder()
one_hot_encoder = OneHotEncoder(categories='auto')
input_features = []
sequences=df.sequences
for sequence in sequences:
     integer_encoded = integer_encoder.fit_transform(list(sequence))
     integer_encoded = np.array(integer_encoded).reshape(-1, 1)
     one_hot_encoded = one_hot_encoder.fit_transform(integer_encoded)
     input_features.append(one_hot_encoded.toarray())

np.set_printoptions(threshold=40)
input_features = np.stack(input_features)
print("Example sequence\n-----------------------")
print('DNA Sequence #1:\n',sequences[0][:10],'...',sequences[0][-10:])
print('One hot encoding of Sequence #1:\n',input_features[0].T)

labels=df.label
one_hot_encoder = OneHotEncoder(categories='auto')
labels = np.array(labels).reshape(-1, 1)
input_labels = one_hot_encoder.fit_transform(labels).toarray()
print('Labels:\n',labels.T)
print('One-hot encoded labels:\n',input_labels.T)

from sklearn.model_selection import train_test_split
train_features, test_features, train_labels, test_labels = train_test_split(input_features, input_labels, test_size=0.25, random_state=42)

def get_batch(x_data, y_data, batch_size):
     idxs = np.random.randint(0, len(y_data), batch_size)
     return x_data[idxs,:,:], y_data[idxs]

epochs = 10
batch_size = 100

x_train = train_features
x_test = train_labels
x_train = x_train / 255.0
x_test = x_test / 255.0
x_test = tf.Variable(x_test)

W1 = tf.Variable(tf.random.normal([784, 300], stddev=0.03), name='W1')
b1 = tf.Variable(tf.random.normal([300]), name='b1')
W2 = tf.Variable(tf.random.normal([300, 10], stddev=0.03), name='W2')
b2 = tf.Variable(tf.random.normal([10]), name='b2')

def nn_model(x_input, W1, b1, W2, b2):
     x_input = tf.reshape(x_input, (x_input.shape[0], -1))
     x = tf.add(tf.matmul(tf.cast(x_input, tf.float32), W1), b1) #ERROR FROM HERE?
     x = tf.nn.relu(x)
     logits = tf.add(tf.matmul(x, W2), b2)
     return logits

def loss_fn(logits, labels):
     cross_entropy = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(labels=labels,logits=logits))
     return cross_entropy

optimizer = tf.keras.optimizers.Adam()
total_batch = int(len(y_train) / batch_size)

for epoch in range(epochs):
    avg_loss = 0
    for i in range(total_batch):
        batch_x, batch_y = get_batch(x_train, y_train, batch_size=batch_size)
        # create tensors
        batch_x = tf.Variable(batch_x)
        batch_y = tf.Variable(batch_y)
        # create a one hot vector
        batch_y = tf.one_hot(batch_y, 10)
        with tf.GradientTape() as tape:
            logits = nn_model(batch_x, W1, b1, W2, b2)
            loss = loss_fn(logits, batch_y)
        gradients = tape.gradient(loss, [W1, b1, W2, b2])
        optimizer.apply_gradients(zip(gradients, [W1, b1, W2, b2]))
        avg_loss += loss / total_batch
    test_logits = nn_model(x_test, W1, b1, W2, b2) ##ERROR LOG APPEARS HERE
    max_idxs = tf.argmax(test_logits, axis=1)
    test_acc = np.sum(max_idxs.numpy() == y_test) / len(y_test)
    print(f"Epoch: {epoch + 1}, loss={avg_loss:.3f}, test set accuracy = {test_acc*100:.3f}%")
print("\nTraining complete!")

Traceback (most recent call last): File "", line 16, in File "", line
4, in nn_model File
"/usr/local/lib/python3.8/dist-packages/tensorflow/python/util/traceback_utils.py",
line 153, in error_handler raise e.with_traceback(filtered_tb) from
None File
"/usr/local/lib/python3.8/dist-packages/tensorflow/python/framework/ops.py",
line 7107, in raise_from_not_ok_status raise
core._status_to_exception(e) from None # pylint:
disable=protected-access
tensorflow.python.framework.errors_impl.InvalidArgumentError: Matrix
size-incompatible: In[0]: [2,1], In[1]: [784,300] [Op:MatMul]

I would check and make sure that the dimensions of your matrix multiplication make sense. Specifically, figure out if you can reproduce the error by just multiplying matrices of the size you're using, since as you point out that seems to be where the error comes from.

The sci-kit learn looks fine, for a proof of concept. Furthermore, the deep learning looks good, e.g. Adams, loss etc .. all good. Of course it at ain't RNN or CNN but its proof of concept.
However, the error you encountered suggests it is about resizing the numpy array and that is better answered on Stackoverflow in the absence of an answer here. Dealing with pandas is one thing, munging a numpy array is pretty hardcore. Its a single line of code, but you've got to know it exactly.
An example is where I d focus,
tf.reshape(x_input, (x_input.shape[0], -1))

Before or whether you need to post there please check:

Where is import tensorflow as tf  ?
Using keras without e.g. from tensorflow.keras import optimizers (well okay you can do that given point 1 but generally the keras method is specifically imported)

If I haven't misread the code it will fall over, as soon as it gets past sci-kit learn.  The advantage of placing all imports at the top of the code is that its easy to read whats going on. There is another issue, but if a big deal Stackoverflow (or someone else here) will flag it.

Just to we explain point 2, I suspect the tensor flow on question is Tensorflow 2 if this is true then  Keras is built in and thus tf.keras is fine

Thank you for the answer. import tensorflow as tf ;  tf.reshape(x_input, (x_input.shape[0], -1)) does not help to solve the error...

Okay fine, just my guess about a point in the code that could generate the error.


text:
Getting error in running GATK
I am trying to making a fasta doc using GATK tools but I don't know what does this error say
Server2:/data2/RNASeq/Angel/gatk$ ./gatk CreateSequenceDictionary -R reference.fa
    Using GATK jar /data2/RNASeq/Angel/gatk/gatk-package-4.2.4.1-local.jar
    Running:
        java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /data2/RNASeq/Angel/gatk/gatk-package-4.2.4.1-local.jar CreateSequenceDictionary -R reference.fa
    INFO    2022-01-28 15:04:19 CreateSequenceDictionary    Output dictionary will be written in reference.dict
    15:04:19.401 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/data2/RNASeq/Angel/gatk/gatk-package-4.2.4.1-local.jar!/com/intel/gkl/native/libgkl_compression.so
    [Fri Jan 28 15:04:19 GMT 2022] CreateSequenceDictionary --REFERENCE reference.fa --TRUNCATE_NAMES_AT_WHITESPACE true --NUM_SEQUENCES 2147483647 --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_secrets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false
    Jan 28, 2022 3:04:19 PM shaded.cloud_nio.com.google.auth.oauth2.ComputeEngineCredentials runningOnComputeEngine
    INFO: Failed to detect whether we are running on Google Compute Engine.
    [Fri Jan 28 15:04:19 GMT 2022] Executing as Server2 on Linux 4.15.0-166-generic amd64; OpenJDK 64-Bit Server VM 11.0.13+8-Ubuntu-0ubuntu1.18.04; Deflater: Intel; Inflater: Intel; Provider GCS is available; Picard version: Version:4.2.4.1
    [Fri Jan 28 15:04:19 GMT 2022] picard.sam.CreateSequenceDictionary done. Elapsed time: 0.00 minutes.
    Runtime.totalMemory()=2155872256
    To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
    picard.PicardException: /data2/RNASeq/Angel/gatk/reference.dict already exists.  Delete this file and try again, or specify a different output file.
        at picard.sam.CreateSequenceDictionary.doWork(CreateSequenceDictionary.java:220)
        at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:308)
        at org.broadinstitute.hellbender.cmdline.PicardCommandLineProgramExecutor.instanceMain(PicardCommandLineProgramExecutor.java:37)
        at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160)
        at org.broadinstitute.hellbender.Main.mainEntry(Main.java:203)
        at org.broadinstitute.hellbender.Main.main(Main.java:289)
    Server2:/data2/RNASeq/Angel/gatk$

picard.PicardException: /data2/RNASeq/Angel/gatk/reference.dict already exists.  Delete this file and try again, or specify a different output file.

I think that requires no further explanation.


text:
(scRNA-seq) I have a dataset with 4 conditions (experiments) should I be anchoring aka integrating data for each comparison when looking for DEGs?
I have an experiment with 4 conditions in it two of them are controls.  When I am correcting for batch effects (either with PCA or CCA) should I restrict my comparisons between particular groups?  So first I should do control 1 with control 2 then hopefully there will be no DEGs between them.  Next I could do Control 2 vs Condition 1 then condition 1 vs condition2 and control 2 vs condtion 2.  Sorry if this is confusing but I am trying to avoid going into the experimental details in the interest of keeping this post brief.
As a sort of follow up question can anyone suggest the "standard" or best seurat method for correcting batch effects.  I've found some tutorials that recommend using CCA and others suggest PCA.

Batch correction/data integration is a tricky subject for scRNA-seq, see a few recent papers for benchmarks: https://www.nature.com/articles/s41592-021-01336-8 and https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1850-9
Batch correction methods may not even impact differential expression. For example, if you correct the PCs, via a popular method like harmony, then this will impact the t-SNE/UMAP plots but will not impact the actual expression values in the counts matrix used to detect DEGs.
For detecting DEGs, the most common methods are simple stats tests like Wilcoxan Rank test that largely ignore batch effects. However, these are likely not the best approaches. I suggest this recent paper (https://www.nature.com/articles/s41467-021-25960-2). Here they state that the best methods for detecting DEGs consider biological replicates, like is done with bulk RNA-seq. By considering replicates, this is effectively like batch correction

I have a total of 16 replicates, 4 per condition or experiment.  I don't understand how I am supposed to correct across those replicates.


text:
Are there different PAM matrices for different protein families? - How is the inferred common ancestor known to start with?
I am confused about the sequences one starts with to generate the PAM matrices.
The sequences to start with are supposed to diverge by no more than 1% (in sequence identity) wrt a common "inferred" ancestor. - if I understood correctly.
How is the common ancestor sequence inferred without having the means to do any sequence alignment to start with?
Does this mean that a PAM matrix is specific for an inferred ancestor, i.e. a protein family? - I got the impression there was only one in practice.
Thanks in advance,
Gordon

Ideally, but I wouldn't worry because the chance of meeting the criteria in practice are slim. Most importantly you need a tree which is 'independently' derived from the protein sequence family in question. The PAM matrix is absolutely dependent on this. Generating a "true" tree without reference to the sequences you are studying isn't trivial, so it's got to generated with a different data set and the question is whether that data set gives a true tree, if you're stuck the DNA sequence of is a rough approximation (not ideal but you could justify it). If you can adequately do this, I wouldn't worry about the remaining assumptions because it's simply an optimisation which can't proven.
The PAM matrix generation will estimate the sequence of a common ancestor, based on a given tree. What this criteria is saying is that if you want a 'perfect matrix' the sequence reconstruction of the ancestral node needs 'x' because that is pivotal to calculating the matrix. If that's a bit out your matrix is a bit by a margin of error that isn't quantified (maybe it has been, see below)
Whether you can get get something close to the 'hub' sequence within your data set is one question and whether that sequence was included in the first place is all trial and error. I would assume the criteria are ascertained has been set by simulation studies (its fairly easy to do if you've got the time), which only affects the margin of confidence within the confines of that given simulation. PAM matrix generation has been around since the dawn of time and founding principle is an independent tree, everything  else you can calculate and arm-wave, ie admit the calculation may not be so precise. In other it doesn't stop you doing it and presenting the results, that is perfectly valid, as long as simulated criteria are understood. The problem would be is if you ignore a past study and present your results 'without fault'. The biggest problem would be if you generated a tree and used the same data.to generate the matrix, because that is whatever matrix you used to generate the tree and that's a no no.

Alright, I see - so in practice, there is only one PAM1 matrix generated from all known protein sequences. (so to speak in the limit... I assume it converges to that)  ---
And when performing any pairwise alignment, I would compare my sequences to each other ensuring that no one pair has more than 1% divergence in identity - instead of comparing to that inferred common ancestor?

Whilst you have the gist, without supervision in coalescent theory this could be tricky. My real concern is the criteria are very theoretical and seldom met in real world data to be meaningful or practical. If If I was in your position I would justify my why I didn't meet the criteria.and move on. The tricky bit is our knowledge of theoretical basis of molecular phylogeny isn't identical and I don't really know your precise question or why it's been set. Given that I'd proceed with caution here, it could be a 'can of worms'


text:
No bootstrap value on single-copy gene tree created by OrthoFinder
I'm running analysis with OrthoFinder and it produces a single-copy gene tree. When I visualize the tree with iTOL, there's no bootstrap value on the branches or nodes. There's a published paper in BMC Genomics that also didn't show the bootstrap value using OrthoFinder.
https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-019-6065-7#Tab1

Does anybody know why there's no bootstrap values present on the phylogeny?
How can this be corrected?

Thanks

Please check the file
 `$PWD/SpeciesTree_unrooted.txt" `

You should find the bootstrap support* values per branch there. If not make sure is the latest version is installed from https://github.com/davidemms/OrthoFinder
The papers are: genomebiology.biomedcentral.com/articles/10.1186/… and genomebiology.biomedcentral.com/articles/10.1186/… Genome Biology 16:157 ...
It officially offers tree stats:
"OrthoFinder also provides comprehensive statistics for comparative genomic analyses"
The documentation is here

The question was raised how to increase replicates by one, two or three logs. This is not really orthodox technique. The way the issue should be first addressed is the user performs replicates manually, e.g. 10 new runs, and then looks across the consistency of bootstrap values for each manual replicate. If all values of interest between manual replicates are >80% or <70% it really doesn't matter. However if a key branch ranges from 70 - 75% that is strictly unresolved and would generally be resolved using an alternative test. If the variation is >= 10% that would need addressing. Generally variation between replicates is very low =< 2%

note, the original answer just contained the word 'support' and has now been more correctly described as "bootstrap support". Thanks to the OP for pointing this out in the comments below.

I see, thanks! If I wanna know the bootstrap number, where could I see it? Can I customize the bootstrap number to something like 10000?

The file I mentioned contains bootstrap values. I have modified the post to provide the link to the documentation where it describes this. It is under the header 'species tree'. Increasing the bootstrap replicaties is has been answered in the modified post.


text:
Bio.codonalign.codonseq module (cal_dn_ds) for SARS-COV-2
I have calculated the dN,dS ratio for a SARS-COV-2 data set ('NC_045512.2' (Wuhan strain) and 'LC666924.1'), and results in a zero value (for most of the genes).
Is it possible and what is the interpretation, for example negative selection?

Yes that is the expected result, because there are no amino acid mutations in the data sets you are using. Thus the numerator is 0 and anything divided by zero is always zero, even if every codon had a synonymous mutation.
The most precise interpretation is simply there are no amino acid mutations. It doesn't matter how you reword that, all divergence is resulting from synonymous mutations.
The only thing to watch out for is a rounding error, i.e. a single amino acid change is dropping off detection because dN/dS is incorrectly rounded.

A post was raised (now deleted) about what does that say about selection. dN/dS theory states that if it is significantly >1 thats positive selection .. etc. However, a truer representation of the data is there is no amino acid mutation, rather than their is a significant excess of synonymous to non-synonymous mutations. A zero value makes things tricky, and is exemplified by two hypothetical data sets one where every codon had one synoymous mutation (strong negative selection) compared against one that had a single synonymous mutation (little mutation) and both data sets had zero amino acid mutations - both data sets would have a dN/dS of zero but that is not a reflection of negative selection.

To answer the point below a single mutation is not a great basis draw statistical significance and the dN/dS ratio is not a good test if either the numerator or denominator are zero.
Most data is under negative selection in strict accordance with dN/dS.
Neutral selection is extremely rare.
Summary Sometimes the statistical analysis of a data set does not say very much ... and sometimes we have to accept that as an answer.

So if that had a little mutation, Does that mean it's a neutral selection?
And if the sequences (codons and amino) are exactly the same? We can say something about the interpretation?
See examples below-
ORF8
Seq1:
MKFLVFLGIITTVAAFHQECSLQSCTQHQPYVVDDPCPIHFYSKWYIRVGARKSAPLIELCVDEAGSKSPIQYIDIGNYTVSCLPFTINCQEPKLGSLVVRCSFYEDFLEYHDVRVVLDFI
Seq2:
MKFLVFLGIITTVAAFHQECSLQSCTQHQPYVVDDPCPIHFYSKWYIRVGARKSAPLIELCVDEAGSKSPIQYIDIGNYTVSCLPFTINCQEPKLGSLVVRCSFYEDFLEYHDVRVVL--I
ORF10
Seq1
MGYINVFAFPFTIYSLLLCRMNSRNYIAQVDVVNFNLT
Seq2
MGYINVFAFPFTIYSLLLCRMNSRNYIAQVDVVNFNLT
E
(The same)
Seq1
MYSFVSEETGTLIVNSVLLFLAFVVFLLVTLAILTALRLCAYCCNIVNVSLVKPSFYVYSRVKNLNSSRVPDLLV
Seq2
MYSFVSEETGTLIVNSVLLFLAFVVFLLVTLAILTALRLCAYCCNIVNVSLVKPSFYVYSRVKNLNSSRVPDLLV

I've answered above

@user18071336 If you have a new question do not ask a question as an answer, post a new question, otherwise you will get downvoted and get banned by the system.


text:
Error in running Mutect2
I am trying to run Mutect2 on a .bam file but I get an error by Googling can not be tackled
Have you seen this error?
fi1d18@RBGO-Server2:/data2/RNASeq/Angel/gatk$ ./gatk Mutect2 -R reference.fa -I sample.bam --germline-resource af-only-gnomad.hg38.vcf.gz --panel-of-normals pon.vcf.gz -O single_sample.vcf.gz
Using GATK jar /data2/RNASeq/Angel/gatk/gatk-package-4.2.4.1-local.jar
Running:
    java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /data2/RNASeq/Angel/gatk/gatk-package-4.2.4.1-local.jar Mutect2 -R reference.fa -I sample.bam --germline-resource af-only-gnomad.hg38.vcf.gz --panel-of-normals pon.vcf.gz -O single_sample.vcf.gz
13:16:43.616 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/data2/RNASeq/Angel/gatk/gatk-package-4.2.4.1-local.jar!/com/intel/gkl/native/libgkl_compression.so
Jan 31, 2022 1:16:43 PM shaded.cloud_nio.com.google.auth.oauth2.ComputeEngineCredentials runningOnComputeEngine
INFO: Failed to detect whether we are running on Google Compute Engine.
13:16:43.733 INFO  Mutect2 - ------------------------------------------------------------
13:16:43.734 INFO  Mutect2 - The Genome Analysis Toolkit (GATK) v4.2.4.1
13:16:43.734 INFO  Mutect2 - For support and documentation go to https://software.broadinstitute.org/gatk/
13:16:43.734 INFO  Mutect2 - Executing as fi1d18@RBGO-Server2 on Linux v4.15.0-166-generic amd64
13:16:43.734 INFO  Mutect2 - Java runtime: OpenJDK 64-Bit Server VM v11.0.13+8-Ubuntu-0ubuntu1.18.04
13:16:43.734 INFO  Mutect2 - Start Date/Time: 31 January 2022 at 13:16:43 GMT
13:16:43.734 INFO  Mutect2 - ------------------------------------------------------------
13:16:43.734 INFO  Mutect2 - ------------------------------------------------------------
13:16:43.735 INFO  Mutect2 - HTSJDK Version: 2.24.1
13:16:43.735 INFO  Mutect2 - Picard Version: 2.25.4
13:16:43.735 INFO  Mutect2 - Built for Spark Version: 2.4.5
13:16:43.735 INFO  Mutect2 - HTSJDK Defaults.COMPRESSION_LEVEL : 2
13:16:43.735 INFO  Mutect2 - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
13:16:43.735 INFO  Mutect2 - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
13:16:43.735 INFO  Mutect2 - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
13:16:43.735 INFO  Mutect2 - Deflater: IntelDeflater
13:16:43.735 INFO  Mutect2 - Inflater: IntelInflater
13:16:43.735 INFO  Mutect2 - GCS max retries/reopens: 20
13:16:43.735 INFO  Mutect2 - Requester pays: disabled
13:16:43.735 INFO  Mutect2 - Initializing engine
13:16:44.015 INFO  Mutect2 - Shutting down engine
[31 January 2022 at 13:16:44 GMT] org.broadinstitute.hellbender.tools.walkers.mutect.Mutect2 done. Elapsed time: 0.01 minutes.
Runtime.totalMemory()=2155872256
***********************************************************************

A USER ERROR has occurred: Cannot read file:///data2/RNASeq/Angel/gatk/pon.vcf.gz because no suitable codecs found

***********************************************************************
Set the system property GATK_STACKTRACE_ON_USER_EXCEPTION (--java-options '-DGATK_STACKTRACE_ON_USER_EXCEPTION=true') to print the stack trace.

I also get this
fi1d18@RBGO-Server2:/data2/RNASeq/Angel/gatk$ file pon.vcf.gz
pon.vcf.gz: HTML document, ASCII text, with very long lines

Sorry, I think the nature of stackexchange is such that you will need to post a new question to get help with the new problem. Perhaps a better location for this question is biostars? https://www.biostars.org/

Yep, edits are for clarification, questions can't be moving targets (I understand it's frustrating, but changing the question renders the answers obsolete and that just generates a mess - which is the reason we don't allow it here)

Your input file pon.vcf.gz should be a gzip file, but from your file command it looks like it is an HTML file instead. Perhaps you downloaded the VCF file incorrectly?
Here is what I get on a bona fide gzip file in Ubuntu:
$ echo hi | gzip -c > hi.gz
$ file hi.gz
hi.gz: gzip compressed data, from Unix, original size modulo 2^32 3

Thank you, you were right, I downloaded the proper file but I get the same error, please have a look at my main post]

Downloading that vcf file from this link fixes the problem (but generates a new one)
wget https://storage.googleapis.com/gatk-best-practices/somatic-hg38/1000g_pon.hg38.vcf.gz


text:
Is use_dimred deprecated or is this a runTSNE() bug?
NOTE: I posted the same question on Biostars but so far no replies there so I figured I would ask here, too.
So, I was running a pipeline which I developed some time ago and a very weird behaviour was found, which had quite an impact over our analyses.
Indeed, the scater's runTSNE() function was subjected to a change in the parameters names and, in particular, the old use_dimred parameter now has a new name, which is dimred.
Interestingly, I was running the pipeline using the runTSNE() function and passing the dimensionality reduction through use_dimred which is apparently deprecated and obtaining very weird results without being aware of that since no warning and/or error was returned.
A quick example:
# load some data
library("scRNAseq")
library("scater")
sce.416b <- LunSpikeInData(which="416b")

check the data and notice no dimensionality reduction is available:
> sce.416b
class: SingleCellExperiment
dim: 46604 192
metadata(0):
assays(1): counts
rownames(46604): ENSMUSG00000102693 ENSMUSG00000064842 ...
  ENSMUSG00000095742 CBFB-MYH11-mcherry
rowData names(1): Length
colnames(192): SLX-9555.N701_S502.C89V9ANXX.s_1.r_1
  SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 ...
  SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1
  SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1
colData names(9): Source Name cell line ... spike-in addition block
reducedDimNames(0):            <-- HERE IT IS EMPTY
mainExpName: endogenous
altExpNames(2): ERCC SIRV

at this point I simply run the runTSNE() function using a non-existent dimred:
# set seed
set.seed(001)
# add logcounts
sce.416b <- logNormCounts(sce.416b)
# add random colours to cells
colLabels(sce.416b) <- rep(c(1,2,3,4), 48)
# runTSNE, first try
sce.416b <- runTSNE(sce.416b, use_dimred="made_up_dimred", perplexity=30)
# and plot
dev.new(); plotTSNE(sce.416b, colour_by="label")
# runTSNE, second try
sce.416b <- runTSNE(sce.416b, use_dimred="another_made_up_dimred", perplexity=30)
# and plot
dev.new(); plotTSNE(sce.416b, colour_by="label")
# runTSNE, third try
sce.416b <- runTSNE(sce.416b, use_dimred="a_third_made_up_dimred", perplexity=30)
# and plot
dev.new(); plotTSNE(sce.416b, colour_by="label")

it, unfortunately, works. Indeed, I now have three TSNE available that I did plot and that are slightly different with each other.
> sce.416b
class: SingleCellExperiment
dim: 46604 192
metadata(0):
assays(2): counts logcounts
rownames(46604): ENSMUSG00000102693 ENSMUSG00000064842 ...
  ENSMUSG00000095742 CBFB-MYH11-mcherry
rowData names(1): Length
colnames(192): SLX-9555.N701_S502.C89V9ANXX.s_1.r_1
  SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 ...
  SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1
  SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1
colData names(10): Source Name cell line ... block sizeFactor
reducedDimNames(1): TSNE      <-- THERE YOU GO: A NICE TSNE COMPUTED OVER A NON-EXISTENT DIMRED
mainExpName: endogenous
altExpNames(2): ERCC SIRV

What is even more interesting is that if I run the function with the newly renamed parameter, i.e. dimred, I do get an error:
> sce.416b <- runTSNE(sce.416b, dimred="made_up_dimred", perplexity=30)
Error in value[[3L]](cond) :
  invalid subscript 'type' in 'reducedDim(<SingleCellExperiment>, type="character", ...)':
  'made_up_dimred' not in 'reducedDimNames(<SingleCellExperiment>)'

But then, if you have an old pipeline which used the runTSNE function with the use_dimred and you're now running an updated version of scater, I have no idea what you may get.
So, as per my question...is this a bug? Has anyone else noticed this weird behaviour? Am I wrong?

I believe it would be easier to answer the question by dividing it into related but more specific questions. Having said that, even if the answers below would be satisfactory, I would still post this on the package's Github repo, the author being informed about this can help other users of this package.
Question 1: Why does runTSNE() work when I use the deprecated use_dimred argument?
Let's create function and call it:
> my_func <- function(dimred) {
+   if(missing(dimred)) {
+     print("dimred argument is not provided so something else is being done...")
+   } else {
+     print(dimred)
+   }
+ }
> my_func(42)
[1] 42
> my_func(dimred = 42)
[1] 42
> my_func()
[1] "dimred argument is not provided so something else is being done..."
> my_func(use_dimred = 42)
Error in my_func(use_dimred = 42) : unused argument (use_dimred = 42)

As expected, we get an error if we call a function with an argument that does not belong to our function: use_dimread. If we define our function with the special argument ... however, it works:
> my_func <- function(dimred, ...) {
+   if(missing(dimred)) {
+     print("dimred argument is not provided so something else is being done...")
+   } else {
+     print(dimred)
+   }
+ }
> my_func(42)
[1] 42
> my_func(dimred = 42)
[1] 42
> my_func()
[1] "dimred argument is not provided so something else is being done..."
> my_func(use_dimred = 42)
[1] "dimred argument is not provided so something else is being done..."

Basically when you call a function, the function you are calling can make calls to other functions and the ... is basically allowing the user to feed arguments for these let's say "secondary" functions. Apparently R does not care if what you feed as an argument to a function is used or not if the function is defined with the special ... argument.
Question 2: Why is runTSNE() working if the supplied use_dimred argument was not used (as shown above)?
If you take a look at the source code, calculateTSNE(), which runTSNE() runs under the hood, is using PCA if the dimred argument is not provided, its default is NULL anyway. It is been a while since I have used Scater but I think this would mean that the tSNE is being computed on (some number of) PCs.
And as a conclusion, I don't think we can really call this a bug but a rather unexpected behavior because of use of the special ... argument in the function definition.

well, this is definitely interesting. I did not consider the implications of using the .... I will let the developers know for sure.


text:
How to identify different protein domains using HHpred
I am new to HHpred/ analyzing proteins so bare with me. I have been given an uncharacterized protein, whose FASTA sequence I have obtained from Uniprot. I am looking to do the following by using HHpred.
So far, I have done the following:

obtained FASTA sequence and inputted it into HHpred.
Loaded the FASTA sequence on to HHpred.

What I am looking to do on HHpred now:

Identify the different domains/motifs of this uncharacterized protein
Compare full length protein to two separate domains
Determine if the predictions has improved by comparing the full length proteins to the separate domains.

All comments help, I am very desperate for any tips.
Thank you!

Do the automated annotation databases, such as PFam etc., in the Uniprot section Family & Domain, have anything annotated? Are they missing the domain you are after (they generally do), is the boundary wrong (they generally are), do you want more granularity (family too big) etc, etc?


text:
Some definitions about RNA-seq
I want to select a 10x single cell kit
What does 2x   50     75     100    150     250 mean in paired end sequencing?

I am pretty sure that is read length. Both reads in a pair (2x) can be 50, 75, ... 250 bases long. And I am even more sure you will find this information on the webpage of the provider and if not, you should email them.


text:
Definition of "proband" in CNV analysis
I am reading this paper regarding de novo CNV analysis in Schizophrenia:

Accordingly, such CNVs must frequently occur as recurrent de novo mutations. In a sample of 662 schizophrenia proband–parent trios, we found that rare de novo CNV mutations were significantly more frequent in cases (5.1% all cases, 5.5% family history negative) compared with 2.2% among 2623 controls, confirming the involvement of de novo CNVs in the pathogenesis of schizophrenia.

I assumed that proband stands for the sample being tested (child of the parents), however, I couldn't find an explicit definition for it. Is it just the sample being tested, or does it imply additional details?

answer from @maximilian-press, converted from comments:

https://www.genome.gov/genetics-glossary/Proband
https://www.merriam-webster.com/dictionary/proband

A proband is an individual affected with a disorder who is the first subject in a study (as of a genetic character in a family lineage) who brings a genetic disorder to the attention of the medical community.


text:
How do I build a glycosylated protein from PDB files manually?
I have a large protein PDB file (~3500 residues) and a couple of PDB files containing glycan structures I sourced from GLYCAM.
How can I attach the N and O-glycans onto each protein residue with software? In PyMol, I am stuck at loading up the protein PDB file and glycan file together, but I don't know how to put the two elements together. Are there any commands I can use? or tutorial resources I can take a look at? I'm open to trying other platforms.
I've tried glycan modelling webtools (GLYCAM, CHARMM-GUI) but they can't seem to process my large glycoprotein - leading to errors or extended waiting times which also end in failure.

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

You can use ChimeraX.

Open your glycans and your protein pdb files as different models, then move the models of your glycans until they are near the residues you want (in the Toolbar: Right Mouse -> Move Model).
Then use combine to combine one of your glycan model with the protein
Finally Use Tools -> Structure Editing -> Build Structure and its Adjust Bonds menu to add covalent bonds. You can also add or delete atoms this way.

This is a straightforward but complicated and tedious way of doing this, and it will give you unideal structures. There should be an equivalent way to proceed using PyMol. For a more automatized way I can only recommend having a look at PyRosetta. There is a Jupyter Notebook tutorial about adding glycosylation.
There may be other tools out there for doing this, like MoleculeKit, but I'm less aware of them. A nasty solution would be to cut your protein in smaller pieces, glycosylate each part with CHARMM-GUI, then reassemble them.

answer from @coot, converted from comments (please edit this answer to improve it if you can):
Manual model building is a common thing in structural biology. Usually, Coot is used.


text:
How to compare across Microarray datasets for a gene in Mice?
How does one go about comparing expression levels of one gene across multiple microarray datasets pulled from different experiments/papers in mice? In my example, I've downloaded the expression levels for my gene of interest from different microarray datasets for different cell populations of interest (specifically healthy hematopoietic cells). How would I pool the datasets together?
Say I have common cell types from multiple datasets such as B-cells, would I be able to just combine them into one data frame and visualize them as one or would I need to normalize some aspect of that?
I'm really confused because I've had different people tell me different things.

If you are after expression of tissue-specific, but otherwise healthy wild-type mice, Bgee is really a great resource for you. They pull a lot of data, curate it and sort it in a single database that easy to browse. They also have lots of tools for comparative genomics (if you would be wondering about expression of an ortholog in a different species).
If you are just after the methodology of pulling expression data together from multiple sources be because you are interested in pulling data from specific papers, you can check their paper about the database. The pipeline is available on GitHub too.


text:
How to solve correlation problems between different samples in scRNA-seq?
I am trying to align and merge different samples from NCBI.
I end up having correlation problem with these sample. The picture below shows an heatmap of the R² by doing a linear regression between 2 samples' count matrices.

Samples A, B and C are from the same project and sample D from another. A, B and C are sequenced using the 10X Genomics for Single Cell V(D)J + 5' Gene Expression, and Chromium Next GEM Single Cell 3ʹ v3.1 for sample D.
All samples have been recovered from NCBI using fastq-dump — internaly threaded to be more efficient — and aligned using cellranger with the default human genome given by 10x.
QC shows some GC content small shifts along side some huge duplication level — even for scRNA-seq.
Am I missing something important when it comes to merge datasets?
Thank you very much

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Thank you, I am trying to update my comment to get it as clear as possible.

Those correlations don't match your description; A looks very different from the others. Can you confirm that the labels are correct for the graph (e.g. show the code snippet that created it)? I've often found that if create heatmap plots in R with labels defined separately from the matrix, the labels appear in the wrong order.

The picture below shows an heatmap of the R² by doing a linear
regression between 2 samples' count matrices.

In addition to doing a heatmap, I'd do a conventional UMAP/t-SNE plot and see how much seperation of the samples there are across all the different clusters/cell types.

Am I missing something important when it comes to merge datasets?

Certainly there are 'technology batch effects' when comparing 3' and 5' data from the same tissue type. You can try to correct for these when merging the dataset by using a batch correction/integration tool (see reviews 1 and 2), but will likely be imperfect.

QC shows some GC content small shifts along side some huge duplication
level — even for scRNA-seq.

If you're using a tool like fastQC then note that high duplication levels are expected for (single cell) RNAseq from highly expressed genes.

Hi,
I did a umap to better view samples' correlation and it appears that the distance between the samples is proportional to the correlation.
As for the technical batch effect, I am looking into it but I find it strange that the three samples A, B and C, which come from the same dataset and use the same sequencing protocol, are so far apart.
I will test those method.
Thanks you very much for your time


text:
Extract sequence context of high-degree nodes in assembly graphs
I often use metaSPAdes to assemble short reads from human microbiomes. My simplified understanding of short-read de Bruijn graph assemblers is that they fail where ambiguous paths cannot be resolved. While it can be said that these points of failure may be due to certain structural features of the underlying sequences -- long repetitive sequences or highly conserved sequences found at multiple loci -- I'm interested in the functional classification of sequences that are proximal to nodes on assembly graphs. Put another way, what kinds of genes are poorly represented in metagenomes because they are poorly assembled? Naively, I expect mobile elements like transposons and integrons to be in such a list of genes, since they are often flanked by repeated sequences and can be present in many copies across multiple genomic contexts, though I want to take an empirical approach.
I've been dabbling with Bandage, which allows you navigate assembly graphs in a GUI, though I'm having trouble extracting features of the underlying network that Bandage builds from the input (.fastg or .gfa). I'd like to know if there are any tools or workflows that would allow me to, given an assembly graph in GFA format ...

filter nodes by degree,
extract sequences proximal to high-degree nodes, and
filter pulled sequences based on sequencing depth, so as to preclude analyzing node-adjacent sequences in the case where the assembler failed due to low depth.

I'm open to the possibility that I'm taking the wrong approach to this problem, so feel free to answer with a frame challenge.

I don't have a recommendation of specific tools, but we might be able to use the GFA file itself to get at this with just bash commands.
For these examples I use the linked GFA here, which is a tutorial with more useful info.
GFA segment (node) lines start with S, link lines start with L (spec here):
S   3438628 CACAATCTAAGAAATCTGAACTACCTGAAACAGGTGGAGAAGAATCAACAAACAACGGCATGTTGTTCGGCGGATTATTTAGCATTTTAG...
L   3461020 +   3461022 +   55M

Links connect the numeric IDs of segments to each other.
Find high-degree nodes
We can just count the occurrences of node IDs in links, which should give us the degree.
egrep "\bL\b" | cut -f4 assembly_graph_with_scaffolds.gfa  | sort | uniq -c | sort | tail

# prints:

      5 3461076
      5 3461246
      5 3461286
      5 591392
      5 606626
      5 634272
      5 721568
      5 766622
      5 822586
      6 3455660

So we can see that there are a lot of links of segment 3455660. We can confirm this in the GFA:
grep 3455660 assembly_graph_with_scaffolds.gfa
# prints:
S       3455660 CTGATAAACGGCGGGAACAGAACCAACACTACGCGTTGCTCCGATCTCAACAACTTCTAAG   KC:i:2936
L       3455660 +       56698   -       55M
L       3455660 +       3455666 +       55M
L       3455660 +       3456290 -       55M
L       3455660 -       3455652 -       55M
L       3455660 -       3456398 -       55M
P       NODE_6_length_300278_cov_129.624456_1   3448698+,3456292-,3448662+,3448670+,3451896-,3460544-,3460388+,3434104-,3434124-,56698+,3455660-,3455652-,3460374+,3460806-,3460824-,3460806-,3460824-,3460806-,3460530-,3449022-,3454210+,3459744-,3459742-    *
P       NODE_8_length_282610_cov_136.723533_1   3460440-,3434818+,3276488+,3461076-,3461070-,3461062-,3461054-,3461046-,3461038-,3460426-,3461020+,3461022+,3461020+,3460516-,3455644+,3455652+,3455660+,3456290-,3448662+,3460546-,3275518-,3387950+,3456304+,3367628+,523826-,3418048-,3456388+,3460376-,3327686+,3455382-,3460586+,655086-,3275436+,3460670+,3278936+,3460218+,3460346+,3460498+,3459664+,3460496+,3459664+,3459200+,3459424-,3458356+,3455426-,3455420-     *
P       NODE_39_length_23557_cov_173.180283_1   3456398+,3455660+,3455666+,3418048-,3273044-,3274776+,3433484-,3433478-,3433470-,3433462-,3456812-,3460134+,3455826-,3319248+,3454700+,3421874+,3453038+,3459348-,3453046-,3234810-,3460596+,3460980+,3460654+  *

Extract sequences proximal to a node
Having extracted the high-degree node, we can use the output of grep above to find links to other sequences:
grep 3455660 assembly_graph_with_scaffolds.gfa | egrep "\bL"
# prints:
L   3455660 +   56698   -   55M
L   3455660 +   3455666 +   55M
L   3455660 +   3456290 -   55M
L   3455660 -   3455652 -   55M
L   3455660 -   3456398 -   55M

We can then go lookup S lines for those segments:
egrep "\bS\W56698" assembly_graph_with_scaffolds.gfa
# prints:
S   56698   GCATGTCGAGAAAACACCTTAGAAGTTGTTGAGATCGGAGCAACGCGTAGTGTTGGTTCTGTTCCCGCCGTT    KC:i:2067

Coverage filtering
Read count (RC tag) is an optional field for GFA, which does not exist in the linked file.
However, we do have k-mer count (KC), which you could use. Presumably high k-mer count would correlate well to coverage. I don't have a suggestion of a nifty tool that filters on this, I'd use a custom script.
Summary
Probably there exists some command line tool that does this stuff, but I don't know about it. It has to all exist inside assemblers, but it might be pretty hard to pull out.
If I had a little more time, I might write a simple tool that uses knowledge of the GFA structure to perform these kinds of data slicing operations, it seems pretty useful. It really shouldn't be hard to do this stuff!
A related approach, that I believe has a somewhat different focus, can be found here.

Thanks Max! This is exactly the kind of insight I am looking for. I read the GFA specification but was honestly confused by it. This makes everything clear. Excellent answer

In my searching, I've also found RGFATools, which provides a framework for parsing and editing GFA files. They include a method for computing the coverage of a segment given the length of a segment, the k-mer count of the segment, and the average k-mer count of all segments. I think I'll try my hand at implementing this for filtering

@acvill that is interesting- I knew about gfapy (https://academic.oup.com/bioinformatics/article/33/19/3094/3883940) though I have never used it. If you find a better way to do it with that rather than all the kludgey shell stuff, consider posting your own answer.

Edit: I've written a short blog post describing my use of nodeSeqs.sh to investigate the functions of sequences at the nodes of short-read assembly graphs for a set of human oral microbiome samples.
https://albertvill.com/posts/nodeseqs_oral/
TL;DR -- transposases and recombinases are enriched, as expected.

Based on Maximilian's excellent answer, I have created a script to address my own question.
https://github.com/acvill/nodeSeqs
Here are the details on how it works:

Link lines (beginning with L) are extracted from the GFA file, and segment occurrences are counted from both the From and To fields. To ensure that the final sequences capture the overlapping (linked) portions of segments, links are processed as segment-orientation pairs. Segment IDs with at least the specified minimum degree (-d) are written to temp_XXXXX/segments.txt, and the complementary segment lines (beginning with S) are pulled. k-mer counts (KC tag) are used to compute coverage for each segment using the approximation from Gonnella & Kurtz 2016; rewritten:

where KS is the k-mer count for a segment, LS is the length of a segment, and k is the k-mer size used in de Bruijn graph construction (-k). For (meta)SPAdes, graphs are built iteratively using an increasing k-mer size, so k is equal to the largest k-mer size used. Sequences are extracted corresponding to segments whose coverage is at least the specified minimum coverage (-m). Sequences longer than the specified context (-c) are truncated to the context length with respect to the link orientation. Including the -t flag will automatically remove the temporary folder when the program is finished, which holds intermediate files and log.txt.

The code is indeed "shell stuff", though hopefully it is not too "kludgey". ;-)

+1 for putting in the k-mer coverage approximation, which I neglected. Not too kludgey, unlike my answer! I think that you could accept your own answer if you liked, it is a better and more useful approach than mine, and fully answers your question.

@MaximilianPress thanks for the feedback. Since your answer provided the framework to create my answer, I'm happy to keep the accept in place.


text:
Galaxy vs Cancer Genomics Cloud
I want to use cloud computing services to analyze my oncological RNA-seq data.
I came across Galaxy and Cancer Genomics Cloud as two cloud computing platforms which are used by the research community.
Can someone give me the pros and cons of using either Galaxy or Cancer Genomics Cloud.

Good question and I get the idea its about cancer, could you state your specific project/requirements please?

So I am trying to analyze RNA Seq data from Fastq files til differential expression of genes. My P.I. wants me to use either Galaxy or Cancer Genomics Cloud from Seven Bridges Genomics

OP update reveals the data is RNA-seq transcript abundance, that is very different from identifying novel driver mutations. So the stuff below isn't relevant.

Just to point out the Broad Institute has been using Google Cloud Platform (GCP) for years and they must be amongst the leaders in processing next gen sequencing data on cancer oncogenesis. They do have specific GCP pipelines available via Google.
I've been tracking their work for sometime because of potential 'cross-over'. It is kinda impressive, albeit there are a few things they (excessively) over complicate.
Just to throw it out there, and it depends what you mean by 'data', if I was processing 'driver mutations' in NGS there would be no question I'd use GCP under the Broad Institutes pipeline (they've changed its name, but you can just email them). So just to labor (labour) the point it depends whatever you mean by 'data'

Thank you so much for that explanation but my P.I. wants me to analyze Fastq files using either Galaxy or Seven Bridges Genomic's Cancer Genomics Cloud.

Sorry can't help, have updated the question though.


text:
What is the difference between flat and stereo stitch id?
For my thesis, I'm working on side effect prediction in a machine learning framework and I am working on data available on SIDER (http://sideeffects.embl.de/download/). If you open the README file, the first description you'll see is:

meddra_all_se.tsv.gz
1 & 2: STITCH compound ids (flat/stereo, see above)
3: UMLS concept id as it was found on the label
4: MedDRA concept type (LLT = lowest level term, PT = preferred term; in a few cases the term is >neither LLT nor PT)
5: UMLS concept id for MedDRA term
6: side effect name

I can not really understand what is the difference between flat and stereo compound ID; is it a chemical descriptor? Is it a notation issue? Moreover, I thought that the two IDS for the same compound could yield to similar results on PubChem but they don't. For example, in the dataset the STITCH ids for the first compound are CID100000085 and CID000010917; on PubChem the first is an aromatic compound and the second is not, they are widely different.
Thanks to everyone who will take the time to help me.

Hey, welcome to Bioinformatics. I see you asked the same question on another SO! Please try to avoid posting the question on multiple communities.

As for your questions, I may be wrong but I think those terms refer to the chemistry 3D representation in space. I think Stereo should stand for Stereoisomers. You can read more about it here. As for flat, it is how you would represent a compound on a 2D space.
The same compound arranged differently in 3D space can be different molecules and thus yield different results.


text:
Plotting a gene in Seurat
I saw in the extensive Seurat documentation for Dimplot (dimensional reduction plot), here, you can plot a gene by specifying it with group.by = "gene" but this does not work in practice.
> p <- DimPlot(DG.combined, reduction = "umap", group.by = "Fos")
Error: Cannot find 'Fos' in this Seurat object
> p <- DimPlot(DG.combined, reduction = "umap", group.by = "Inhba")
Error: Cannot find 'Inhba' in this Seurat object

Am I mistaken?  Do I need to add a metadata column containing expression of the gene?  How can I access individual cells gene expression data in seurat so I can create a metadata label based on binned expression levels?  Also how would I tell seurat to use a color gradient to label them?
References
https://satijalab.org/seurat/reference/dimplot

As @swbarnes2 has suggested, you can plot gene expression over the top of a dimensional reduction using the FeaturePlot function, using the features parameter for gene names:
p <- FeaturePlot(DG.combined, reduction = "umap", features = c("Fos", "Inhba"))

See the 3k tutorial here:
https://satijalab.org/seurat/articles/pbmc3k_tutorial.html#cb49

I am not totally sure what you are trying to do, but I think you want FeaturePlot, not DimPlot


text:
tblastn not providing full alignment even on exact-matching sequences
I am using tblastn to search a polypeptide seq against transcripts.  One of those transcripts contains the residues for that exact polypeptide once translated, yet BLAST only gets 88% of it rather than it being full length.
Here's the peptide:
>Syntaxin
MTKDRLAALKAAQSDDDDNDDVAVTVDSSGFMEEFFEQVDEIREMIDKIASNVDEVKKKHSAILSAPQTDD
KMKEELEELMSEIKKNANKVRAKLKVIEQNIEQEEHTNKSSADLRIRKTQHATLSRKFVEVMNDYNACQID
YRERCKGRIKRQLAITGKTTTNEELEDMIESGNPAIFTQGIIMETQQANETLADIEARHNDIMKLETSIRD
LHDMFMDMAMLVESQGEMIDRIEYNVEQAVDYIETAKMDTKKAVKYQSKARRKKIMILVCLAILIIILVGV
IGGTLG

And the transcript.
Here's the first-frame translation of the transcript which contains the exact polypeptide above:
frame +1
> - 669
PATRNITIVSIYDYH*NNDLKCVFTE*MEILNRIWILLNSHSP*AQRLHAI**MTKDRLAALKAAQSDDDD
NDDVAVTVDSSGFMEEFFEQVDEIREMIDKIASNVDEVKKKHSAILSAPQTDDKMKEELEELMSEIKKNAN
KVRAKLKVIEQNIEQEEHTNKSSADLRIRKTQHATLSRKFVEVMNDYNACQIDYRERCKGRIKRQLAITGK
TTTNEELEDMIESGNPAIFTQGIIMETQQANETLADIEARHNDIMKLETSIRDLHDMFMDMAMLVESQGEM
IDRIEYNVEQAVDYIETAKMDTKKAVKYQSKARRKKIMILVCLAILIIILVGVIGGTLG***PSSLPLVYS
LIQLFYAACTTAASNIWTHHNRNSYYYNSNSNSNSNCYSNYFSNYYSNRKHYYYYHYYYYYYYYYYYYYYY
YYYYYYDYYYYYYYYYYYYYYYYYCYYYYYYYYYHFHYYYHHYYYHYYYYHYYYHYHHHHHYHYHYHYHHH
HHYHYHYHYHYHYYYYYYYYYYYYYYYYYYYYYYCYYYYYYTFSLVISLYGSSEM*MRLLVLSLYTFYGLF
YPLQQKGSFFF*SLQQTGTGFYNFFSYCQLHNLSIPSPSYFF*EWVFSFLKKHLIEFFCFLILLKFVLSFS
PIISLIISVCRKKISRIFYTIYLKMS**NHYLLF

Yet when I do a local tblastn (using the modern blast+) or use NCBI's pairwise blast tool it doesn't return a full-length match.
Why does it leave off the last 12% of the sequence match?
Here's the alignment:

That's because by default, blast masks low complexity regions. With the online blast2seqs tool, you can disable this by clicking on the "Algorithm options" button:

And then deselect the "Low complexity regions" filter:

This will let you align the entire sequence as expected:
Alignment statistics for match #1
Score   Expect  Method  Identities  Positives   Gaps
Frame
588 bits(1516)  0.0 Compositional matrix adjust.    290/290(100%)   290/290(100%)   0/290(0%)
+1
Query  1     MTKDRLAALKAAQSDDDDNDDVAVTVDSSGFMEEFFEQVDEIREMIDKIASNVDEVKKKH  60
             MTKDRLAALKAAQSDDDDNDDVAVTVDSSGFMEEFFEQVDEIREMIDKIASNVDEVKKKH
Sbjct  154   MTKDRLAALKAAQSDDDDNDDVAVTVDSSGFMEEFFEQVDEIREMIDKIASNVDEVKKKH  333

Query  61    SAILSAPQTDDKMKEELEELMSEIKKNANKVRAKLKVIEQNIEQEEHTNKSSADLRIRKT  120
             SAILSAPQTDDKMKEELEELMSEIKKNANKVRAKLKVIEQNIEQEEHTNKSSADLRIRKT
Sbjct  334   SAILSAPQTDDKMKEELEELMSEIKKNANKVRAKLKVIEQNIEQEEHTNKSSADLRIRKT  513

Query  121   QHATLSRKFVEVMNDYNACQIDYRERCKGRIKRQLAITGKTTTNEELEDMIESGNPAIFT  180
             QHATLSRKFVEVMNDYNACQIDYRERCKGRIKRQLAITGKTTTNEELEDMIESGNPAIFT
Sbjct  514   QHATLSRKFVEVMNDYNACQIDYRERCKGRIKRQLAITGKTTTNEELEDMIESGNPAIFT  693

Query  181   QGIIMETQQANETLADIEARHNDIMKLETSIRDLHDMFMDMAMLVESQGEMIDRIEYNVE  240
             QGIIMETQQANETLADIEARHNDIMKLETSIRDLHDMFMDMAMLVESQGEMIDRIEYNVE
Sbjct  694   QGIIMETQQANETLADIEARHNDIMKLETSIRDLHDMFMDMAMLVESQGEMIDRIEYNVE  873

Query  241   QAVDYIETAKMDTKKAVKYQSKARRKKIMILVCLAILIIILVGVIGGTLG  290
             QAVDYIETAKMDTKKAVKYQSKARRKKIMILVCLAILIIILVGVIGGTLG
Sbjct  874   QAVDYIETAKMDTKKAVKYQSKARRKKIMILVCLAILIIILVGVIGGTLG  1023

This is the part of your transcript sequence that was masked:
acgctaggctgaccctcctccctccctctagtctactccttaatacaattgttctacgccgcgtgcacc

I am guessing it is the repeated ccct nucleotides that resulted in the sequence being masked, but I am not sure.


text:
Nextflow: dynamically setting resource directives for a process
I'm trying to test variations of CPU and memory for a specific process (samtools sort) to measure computational efficiency (measured only by time). To do this, I'd like to set the CPU and memory for a Nextflow process dynamically.
Ultimate question: How can I dynamically set CPU and memory for many instances of a single Nextflow process definition?
Attempt 1
The most obvious solution was to submit distinct Nextflow workflows for each combination of CPU and memory by specifying the CPU and memory on the command line (i.e., when calling nextflow run ...), where my Nextflow workflow has global parameters params.samtools_total_threads and params.samtools_mem_per_thread.
This didn't work perfectly because there appeared to be 'collisions' between parent Nextflow instances trying to lock the history file. I got the following error:
Can't lock file: /path/to/project/.nextflow/history -- Nextflow needs to run in a file system that supports file locks

I could try delaying submissions, but I suspect that will still result in collisions.
Here's the SLURM script that I'm using:
TOTAL_THREADS=$1
MEM_PER_THREAD=$2

nextflow run 02-SAMTOOLS_SORT.nf \
    -with-report "reports/report.${TOTAL_THREADS}_threads.${MEM_PER_THREAD}_GB_per_thread.html" \
    -with-trace "traces/trace.${TOTAL_THREADS}_threads.${MEM_PER_THREAD}_GB_per_thread.txt" \
    --samtools_total_threads $TOTAL_THREADS \
    --samtools_mem_per_thread $MEM_PER_THREAD

My workflow had the following global parameters:
/*
 * Make this pipeline a nextflow 2 implementation
 */
nextflow.enable.dsl=2

...

params.samtools_total_threads = '4'
params.samtools_mem_per_thread = '8'

...

Attempt 2
After Attempt 1 failed, I figured the most 'Nextflow-like way' would be to generate tuples for each combination of CPU and memory that I want to test. Unfortunately, it doesn't appear that it's possible to set CPU and memory directives using parameters passed into the process, but only with global variables.
I get the following No such variable error with the below code (note: most of this code is courtesy of @Steve):
N E X T F L O W  ~  version 21.04.3
Launching `02-SAMTOOLS_SORT.nf` [desperate_aryabhata] - revision: fb9e79055a
No such variable: total_threads

 -- Check script '02-SAMTOOLS_SORT.nf' at line: 33 or see '.nextflow.log' file for more details

Process code
process samtools_coordinate_sort_proc {

    tag { "${bam.baseName}" }

    executor='slurm'
    queue = 'normal'
    cpus = total_threads
    memory = total_mem.GB
    clusterOptions = "--time 5:00:00 --account <account>"

    /* delete files upon completion (I think) */
    cleanup = true

    input:
    path(bam)
    val(total_threads)
    val(total_mem)

    output:
    path("${bam.baseName}.csorted.bam")

    script:

    def additional_threads = task.cpus - 1

    def avail_mem = task.memory ? task.memory.toGiga().intdiv(task.cpus) : 0
    println "avail_mem: $avail_mem"
    def mem_per_thread = avail_mem ? "${avail_mem}G" : ''
    println "mem_per_thread: $mem_per_thread"

    """
    samtools sort \\
        -@ "${additional_threads}" \\
    -m ${mem_per_thread} \\
        -o "${bam.baseName}.csorted.bam" \\
    -T "${bam.baseName}.csorted" \\
        --write-index \\
        "${bam}"
    """
}

Possible attempt 3
Of course, it would be super easy to do this with a couple of bash scripts, but I'm trying to do it in Nextflow for two reasons:

To learn Nextflow better
Because the traces provided by Nextflow are super convenient (though I could get this information from SLURM, too).

It's possible to specify dynamic computing resources using a closure. For example, this sets the cpus and memory directives using the values supplied via the input declaration:
nextflow.enable.dsl=2

process test_samtools {

    echo true

    cpus { cpus }
    memory { memory }

    input:
    tuple path(bam), val(cpus), val(memory)

    script:
    def additional_threads = task.cpus - 1

    def avail_mem = task.memory ? task.memory.toGiga().intdiv(task.cpus) : 0
    def mem_per_thread = avail_mem ? "-m ${avail_mem}G" : ''

    """
    echo samtools sort \\
        -@ "${additional_threads}" \\
        ${mem_per_thread} \\
        -o "${bam.baseName}.csorted.bam" \\
        -T "${bam.baseName}.csorted" \\
        --write-index \\
        "${bam}"
    """
}

workflow {

    cpu_ch = Channel.of( 1, 2, 3, 4, 5)
    memory_ch = Channel.of( 1.GB, 2.GB, 4.GB )

    Channel.fromPath( './data/*.bam' )
        | combine( cpu_ch )
        | combine( memory_ch )
        | test_samtools
}


text:
plot the relation between RMSD and free energy
how can i plot the relation between rmsd and free energy like the image i attached as i need to find the most stable conformation structure for my cycle peptide after finishing REMD simulation

I look at that graph and I note the high ∆G at 0 Ångström and wonder if the reference pose for the RMSD was not energy minimised correctly... so maybe not the best reference. http://prody.csb.pitt.edu/ is a helpful toolkit for this kind of analysis. Is Python an option and what is your reference? (If any user knowledgeable about peptides sees this please step in as I doubt I can help —NMR works by magic, right?)

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Thanks for the comment .. i'm trying to use the tool but couldn't .. in which part please i can plot RMSD with energy surface ?

You would need to calculate the free energy along a trajectory and the rmsd along that trajectory and then cross-correlate/eliminate variables.  There may be a utility which does this but I am not ware of it.
https://manual.gromacs.org/current/user-guide/cmdline.html
There are two utilities of interest for doing just that.  See gmx_energy and gmx-rmsd.  Then you can plot one against the other by elimination of the time/space variables. Check out other calculations you can do with GROMACS:
https://manual.gromacs.org/current/user-guide/cmdline.html.


text:
Querying The Cancer Genome Atlas (TCGA) for gene expression
The UCI gene expression cancer RNA-Seq dataset (https://archive.ics.uci.edu/ml/datasets/gene+expression+cancer+RNA-Seq) consists of 5 types of cancers (BRCA, KIRC, LUAD, COAD, PRAD) and gene expression levels of 20,531 genes for 801 samples. The dataset looks like this:

How to query in The Cancer Genome Atlas (TCGA) to create an updated version of this dataset?
I've done the following query.

Is the given query and the UCI dataset equivalent? If not, how to query properly?
P.S.: I'm a CS student; do not have any biology background.

Not quite what I was hoping for the query rather than the raw link ... Anyway I guess I wasn't clear

To start with, your query above is selecting for breast cancer samples (TCGA-BRCA) only. The original dataset seems to contain 4 other cancer types.
The link you provided has a link to https://www.synapse.org/#!Synapse:syn4301332 and by looking at filenames of the 14 files there, I can tell that there are different types of data: miRNA expression, gene expression, methylation, ... And your query is limited to gene expression only. I am not even sure if you can use just a single query to fetch all of the 14 data files, you might need to use multiple queries.
If you are familiar with R, I would recommend TCGABiolinks, it will abstract away a lot of otherwise annoying manual steps. Apparently you can also use Python retrieve TCGA data:https://docs.gdc.cancer.gov/API/Users_Guide/Python_Examples/

I've added an image of the UCI dataset. I just need to perform a query that returns similar kind of gene expression level tabular data.
As the query only returns TCGA-BRCA data only, will replicating the query for other four types create similar (biologically equivalent and meaningful) kind of dataset?

Yes, you can change the cancer type and re-use the query but my main concern is the type of data that this query would return. Your query fetches count data via the HTSeq - Counts filter but in the posted image, there are floats instead of integers (gene counts), this can be normalized data. If I remember correctly, HTSeq-Counts should be returning actual counts, so integers.


text:
Programmatic way of accessing the FASTA sequence similarity tool (or similar) in Python
I am looking for a tool that performs a sequence similarity algorithm for proteins. More specifically, I am looking for something that would be usable in Python (or anything else usable in the command line actually) which has a similar output to the FASTA tool, i.e:

For my use case, the 'Source' column is not super relevant, so it does not matter if it's not included in the output. What's most important is the DB:ID and the match statistics.
Have you used anything similar? Do you have any suggestions regarding such a tool?
Thanks in advance!

Hmmm ... two suggestions, you can do both via Python, Blast or a tree matrix. Its an interesting question but difficult to know what you're looking for

Will you be comparing specific, known sequences, or is part of the question finding similar sequences to your query in a large database of non-related sequences? Basicaly, do you want a sequence aligner (mafft, t_coffee, clustal etc) or tools like the various flavors of blast (psi-blast, p-blast)? Or are you looking for proteins with shared domains which means tools like HMMER or psi-blast again? Pretty much everything in bioinformatics can be done from the command line, so you need to explain exactly what question you're trying to answer for us to suggest the better options.

Apparently there is a REST API for FASTA. Looks like there are few dependencies but other than that should be straightforward to use.


text:
Understanding fast5 format
I would like to reproduce the results on figure 1 from THIS publication.
I downloaded the data and explore it with h5ls and h5dump based on HDF5.
I would like to understand:

h5ls -r <_fast5file_> | grep Dataset:

/Analyses/Basecall_2D_000/BaseCalled_2D/Alignment Dataset {6925}
/Analyses/Basecall_2D_000/BaseCalled_2D/Fastq Dataset {SCALAR}
/Analyses/Basecall_2D_000/BaseCalled_complement/Events Dataset {5576}
/Analyses/Basecall_2D_000/BaseCalled_complement/Fastq Dataset {SCALAR}
/Analyses/Basecall_2D_000/BaseCalled_complement/Model Dataset {4096}
/Analyses/Basecall_2D_000/BaseCalled_template/Events Dataset {5793}
/Analyses/Basecall_2D_000/BaseCalled_template/Fastq Dataset {SCALAR}
/Analyses/Basecall_2D_000/BaseCalled_template/Model Dataset {4096}
/Analyses/Basecall_2D_000/HairpinAlign/Alignment Dataset {5107}
/Analyses/Basecall_2D_000/Log Dataset {SCALAR}
/Analyses/Calibration_Strand_000/Log Dataset {SCALAR}
/Analyses/EventDetection_000/Reads/Read_0/Events Dataset {11465}

Alignment? Then complement/Fastq is not the complement of template/Fastq (analyzed with h5dump). What is it? Also, events have sometimes 4

Some other datasets that were acquired somewhere else have the raw data and fastq but there is no associated information about kmer and mean(current)/std/length/etc... How do I generate events from raw data?

I went through all documentations and tutorials that I could found in the web but couldn't find anything simple and complete. It feels as if it was a secret.
Example fast5file would be ecoli_er2925.MSssI.timp.100215.fast5/pass/kelvin_100215_ecolisssI_2515_1_ch195_file0_strand.fast5 from PRJEB13021

FAST5 is a proprietary format developed by Oxford Nanopore Technologies, and there's not much great documentation online. All FAST5 files will have the Raw/ field, which contains the original measured Signal. Additional Analyses/ fields can be added by tools such as basecallers which convert signal->ACGT (Guppy) or signal-level aligners (Tombo).
In particular, /Analyses/Basecall_2D_000/ contains the results of 2D basecalling, added by Guppy when run with the --post-out flag. I think Tombo adds the event-level annotations you mention (k-mer mean/std/length). In 2D basecalling, both the physical forward and reverse strands for the same stretch of DNA are measured and basecalled. The basecalls will contain errors, so these will not be exact reverse complements of one another.
Hope this answers your questions, I have a slightly more in-depth explanation on my blog: https://timd.one/blog/genomics/file_formats.php

FWIW, this particular FAST5 was generated in 2017, before Guppy was available

I'll give a different perspective. Oxford Nanopore FAST5 is NGS within a HDF5 database. A clear description is here.
HDF5 is a format used for rapid I/O of large data sets.
Oxford Nano appear to have written a bit boilerplate around HDF5 and called it something that sounds like 'fasta'.
A cool part of HDF5 (not FAST5) is  HDF5 View, which can be downloaded here. Its cool because you can visualise the file structure of HDF5 db via a GUI. I very much suspect if you use HDF5 View on your FAST5 database you will be able to navigate (see not change) the entire file structure.
To answer your question directly h5ls is listing the HDF5 database. To construct a HDF5 db you assign the data under a series of 'keys' which forms the directory structure. However, you perhaps might not need Oxford tools, because  Python pandas have direct series of methods to connect, construct and manipulate HDF5.
I raised the HDF5 in a question then answered here on Bioinfo SE here.
Again the reason for HDF5's importance is unlike an SQL connection, data transfer and disconnect, HDF5 shifts data very quickly.
I do question how Oxford Nanopore can lay claim to FAST5 because HDF5 is open source and always has been. If they want to lay claim to some boilerplate, be my guest - I write the stuff all the time ;-)
HDF5 however is a serious project with major development time and appears to still be under active development.


text:
Questions about microbial 16s rRNA phylogenetic, ANI Taxonomy, and GGDC
I have a few questions about these three methods.

We commonly use 16s rRNA to identify species and construct a phylogenetic tree. If there's a new isolate from a new and undiscovered ecosystem, should we confirm this new isolate identity with ANI taxonomy/ genome to genome distance calculation (GGDC)?

If that isolate has a high similarity with a certain species in the 16s rRNA analysis result, should we proceed to analyze further with ANI/GGDC? Or just let it be because we have a robust 16s rRNA analysis? How can we be sure with our 16s rRNA BLAST result?

When should we use 16s rRNA rather than ANI/GGDC and vice versa?
Is it always necessary to confirm 16s rRNA analysis with ANI/GGDC? (Assuming there's a lot of closely-related genome data available in databases)
Should we confirm the ANI result with GGDC and vice versa?

I do understand that you are working with  is an isolate and thus what occurs in metagenomics you feel is over doing it, particularly if 16S is e.g. >99% similarity via Blast.
Before reading further you might  check my answer history to understand the angle I am approaching this question from. I would see select given methods as the 'gold-standard' other methods are compared against, in addition to single versus multi-gene trees. There is good theoretical criteria for believing this. For example, bacteria are often subject to heavy %GC bias and this makes other methods error prone, because it introduces compositional bias (below).
Summary
Blast alone - thats a 'no' (theory section).  The problem with any one approach is that errors can arise from many different sources stemming from both biology and methodology. Thus a combination of approaches is required. Therefore your maximum approach is a minimum standard.
I personally would go one step further than your question: full maximum likelihood and Bayesian phylogenetics with bootstraps in addition to Blast (minor), ANI and GGDC.
History
From tradition 16S is the major standard, but even then multi-gene phylogenies were required. Thus genome to genome distance calculation (GGDC) with ANI is essential at a minimum. 16S  is back in fashion because its easy to perform a metagenomics study of species population dynamics.
Species delineation has always been problematic in bacteria and the rule of thumb was 5% divergence defined a species boundary.
Advice
I would use a multi-locus appproach under a full phylogenetics analysis, where a species boundary is defined by >5% divergence and bootstraps >80% with Bayes P values > 0.95.
Recombination is essentially ignored in taxonomy providing there isn't over reliance on a single gene tree, but using BLAST, ANI, GGDC as a method data-mining would justify a select sample size. The classification is then nailed with a phylogenetics analysis.
One of your methods (probably ANI) might be based machine learning/deep learning (its not GGDC) and the accuracy of misassigning a species under this method is high. Historically it was really high (80% accuracy), this might have improved to 90%. Its not like image recognition where the error rate is beyond human perception.

Theory A fundamental issue is that recombination in bacteria exists on a spectrum from clonality to panmixia. This makes single gene taxonomy risky, thus at theoretical
level would not believe 16S was a singularly criterion, but it never was. Your bacteria is unknown (to us) so we'd have to assume some level of recombination.
Compositional bias is well known for taxonomic misassignment. It is readily demonstrated in simulation. Thus if two unrelated bacteria demonstrate e.g. 75% GC, they are more likely to artifactually grouped together on e.g. a distance criterion.

Agree. 16S is not reliable for classification, especially at the species level. Bacteria can be quite promiscuous with horizontal gene transfer, so a single gene on its own is not a reliable indicator of taxonomy.

Cool, I did check ANI and its a pairwise aligner and identity algorithm not machine learning. It looks locked behind a webserver so you'd need e.g. Blast to formulate the subset of bacteria to check.


text:
What is the correct order of flooring-normalization-batch correction for microarrays?
This question was also asked on Biostars
I am trying to learn and understand the correct order of data processing steps for microarrays.
I have data which already was analyzed by a researcher using Partek (with an ANOVA model with batch effects). I want to test the results versus a specific gene set in GSEA. However, as I was informed, the genes were not filtered for low expression before ANOVA. So I want to filter genes out with low expression as is recommended for GSEA. However, I wonder whether it is legitimate to filter after DE was already done.
In what order should one apply flooring (setting all values below a threshold to that threshold's value), normalization, and batch correction?

The most important thing to do first is batch correction, and it looks like that was done in this instance.
The ordering in this case seems reasonable. Normalisation is usually applied to a full dataset (it depends on the assumed model that was used for normalisation, and I'm not familiar with how Partek does things).
Filtering after DE should be fine. I'm at least aware that low-count filtering with DESeq2 [cDNASeq rather than microarray] shouldn't influence differential expression results, but it can speed up the calculations slightly:

While it is not necessary to pre-filter low count genes before running the DESeq2 functions, there are two reasons which make pre-filtering useful: by removing rows in which there are very few reads, we reduce the memory size of the dds data object, and we increase the speed of the transformation and testing functions within DESeq2.

https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#pre-filtering

I think that in DESeq, low-count filtering is done prior to normalization so it can influence differential expression slightly.  But that is another topic.


text:
Where to get BED files for regions that contain binding site motifs for specific transcription factors?
Quick question, is there a place where one can download BED files for the binding motifs of certain transcription factors. In my case, I'm looking for a BED file for the genomic regions (enhancers) where TET2 binds to genome-wide. I'm looking to examine the DNA methylation levels in those regions.
Thank You!

The standard answer would be to search JASPAR, to get a curated dataset. However, there does not appear to be anything coming up for TET2 in the database.
If your TET2 gene is what you want specifically, you can find papers that perform ChIP-seq and infer binding regions. I believe that this Supp Table 3 contains information equivalent to a BED file for TET2.