text:
How to statistically test the strength and weakness of a correlated pair?
Subset of df before groupby looks like
Diabetes_Group  DM-1    DM-2    corr
Group1  56  56  1
Group1  56  58  0.611
Group1  56  61  0.448
Group1  56  63  0.3587
Group1  58  56  0.611
Group1  58  58  1
Group1  58  61  0.60113
Group1  58  63  0.534858
Group1  61  56  0.448
Group1  61  58  0.60113
Group1  61  61  1
Group1  61  63  0.622206
Group1  63  56  0.3587
Group1  63  58  0.534858
Group1  63  61  0.622206
Group1  63  63  1
Group2  78  78  1
Group2  78  85  0.622703
Group2  78  94  0.622244
Group2  78  23  0.687826
Group2  85  78  0.622703
Group2  85  85  1
Group2  85  94  0.697313
Group2  85  23  0.612089
Group2  94  78  0.622244
Group2  94  85  0.697313
Group2  94  94  1
Group2  94  23  0.603834
Group2  23  78  0.687826
Group2  23  85  0.612089
Group2  23  94  0.603834
Group2  23  23  1
Group3  27  27  1
Group3  27  80  0.118955
Group3  27  147 0.32038
Group3  27  9   0.335264
Group3  80  27  0.118955
Group3  80  80  1
Group3  80  147 0.430287
Group3  80  9   0.406426
Group3  147 27  0.32038
Group3  147 80  0.430287
Group3  147 147 1
Group3  147 9   0.546452
Group3  9   27  0.335264
Group3  9   80  0.406426
Group3  9   147 0.546452
Group3  9   9   1

What statistical test can be used to do comparative analysis between the different ID's (DM1 or DM2) to check which ID pair has weaker/stronger correlation within each diabetes_group?

Is there any statistical test to categorise the data into three categories:
(a) pairs of DM1 AND DM2 which are highly correlated than others.
(b) pairs that are not highly correlated
(c) For the category 1 pairs (highly correlated), is there any particular ID that outperforms than other ID's among different groups?

This comprehensive analytical outline for the stats and all the stats questions could be captured inside 3 analyses. Implementing this will obviously take a bit of time, but there should not be any hidden caveats (transformations etc ...).
Question 3
The statistical test you want is K-means clustering based on the correlation score where the number of predefined groups is 3. The labels ... you will need to combine the group with the
Label          Value
Group1-56-56   1
...

This is an absolute solution to the data and really cool use of K-means.
Question 3c
An easy solution is:

perform a new K-means with e.g. 10 groups and find the group with highest correlation values.

A better solution is:

hierarchical clustering would identify specific groups and for this sub-question is a good analysis. The input is the correlation matrix, the output is a tree. HOWEVER, interpreting this could be difficult if you're not trained in reading trees.

Question 2
A 2-way ANOVA would be a really good solution to this data for Q2. The test would be the columns were DM1 DM2 and the rows Group1, Group2, Group3. Interpreting it would be fun because you'd also have to grapple with a concept called "interaction" - which could be important here. However, the ability to address a whole raft of questions within one analysis is what 2-way ANOVA would permit. There is a slight issue of unequal group sizes 1, 2 and 3, because it would be better if Group 1,2 and 3 etc... were equal in size. However, the raw data looks to be very amenable to ANOVA. You must not use the correlation index in this test - just the raw values.

K-means clustering is in scipy and scikit-learn
Hierarchical clustering is also in scipy
For ANOVA
import statsmodels.api as sm
from statsmodels.formula.api import ols

There will be easy approaches to 2-way ANOVA e.g. SPSS and I think Excel (needs a plug in). I suspect that Hierarchical clustering will also be in a nice point and click package like SPSS (I don't use them, so I don't precisely know.)

Is it possible for you to explain this in terms of code or any example code

Hi @Megha, you have asked an entire statistical project within a single question, so I answered as best I could.What I propose is to investigate the facets (K-means is a good start) and have a go and ask a specific question if you're stuck. I've updated the post with the location of the specific packages. You don't need to stick with Python, point and click packages are perfectly fine because this is end analysis. Beyond this there is no further process in a data pipeline


text:
How to find specific types of assemblies for specific species using entrez tools?
How to find specific types of assemblies for specific species using entrez tools?
Task: Trying to specifically find transcriptomes and associated cDNA data for a list of speices.
I can use this oneliner esearch -db assembly -query "{species name }" | esummary to get a list of assemblies and associated data available for any given species.
I can use esearch -db assembly -query "{species name }" | xtract -pattern ENTREZ_DIRECT -element Count to find any specific info nested within the XML tree.
Now I am not sure how to find if the assembly is a transcriptome (which is what I am looking for). I have tried to CRTL +f terms like rna, peptide,amino acids to no avail and I guess there is a chance I am missing something while looking it over manually.
For example here is the output for "zea mays".
Does anyone have any clues? Or should I be using a different db or different tool?

I strongly recommend switching to SRA https://www.ncbi.nlm.nih.gov/sra
Whether the following would work ...
esearch -db sra -query "{species name }" | esummary

I'm not at all sure because SRA has its own toolkit. Maybe NCBI e-utilities e-summary efetch would work on the SRA database. However, the weblink certainly works for transcriptomics.
NCBI has shifted some of its database to be housed off-site onto GCP and AWS and I have a feeling these may not be accessible to e-utilities, but I could be wrong.
Just to note, what I envisage the question meaning is a complete transcriptomic study for a given species (or given sample) under given experimental conditions.

Comments. I better clarify, SRA will be raw Illumina output, it is not assembled into exome cDNA. e-utils is basically for assembled cDNA, it might give you the information - I really don't think it will download it though. So to answer your questions, yes the data will be there on SRA, whether this is the format you want is a different question.

I am actually looking for RNA-seq data for the whole exome. I will try to see if I can find transcriptions and cDNA using SRA. Hopefully off-housing has not affected the API.

@Sudoh Exome, that would be assumed .... there's a format issue. Post-updated

source: Stack Exchange

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

sort: votes


text:
Velvet Optimizer automatically changes to hash-length 31
I'm trying to use Velvet Optimizer for a De Novo Assembly; I set my hash-lengths to be between 55 and 69
/projects/path/to/perlscripts/VelvetOptimiser-2.2.6/VelvetOptimiser.pl -s 55 -e 69 -f '-shortPaired -fastq.gz -separate trial2.trimmed.paired.R1.fastq.gz trial2.trimmed.paired.R2.fastq.gz' -f '-short -fastq.gz trial2.trimmed.unpaired.R1.fastq.gz' -f '-short -fastq.gz trial2.trimmed.unpaired.R2.fastq.gz' -x 2

When I check the output file, it returns me with
Starting to check input parameters.
    Start hash value too high.  New start hash value is 31.
    End hash value not in workable range.  New end hash value is 31

Any idea where I could be going wrong in my command?

If a value for MAXKMERLENGTH was not set when compiling velveth, VelvetOptimizer just uses the default $maxhash value of 311. The hash lengths that you specify on the command line must be less than or equal to this value to run the tool successfully. The solution is to re-compile velveth with a larger MAXKMERLENGTH, for example:
make 'MAXKMERLENGTH=69'

Alternatively, re-install using the bioconda package manager:
conda install -c bioconda perl-velvetoptimiser

This package includes velvet, which is built using a much higher maximum k-mer length and is usually sufficient for most purposes: MAXKMERLENGTH=1912.

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


text:
To get background_gradient plot for all data frames in the for loop
I tried the below code to get the background_gradient plot for each diabetes_group. It is giving me only plot for the last group in the plot. How to get background_gradiant for all items (x) in diabetes_group with item name on the top ?
Subset of df before groupby looks like
Diabetes_Group  DM-1    DM-2    corr
Group1  56  56  1
Group1  56  58  0.611
Group1  56  61  0.448
Group1  56  63  0.3587
Group1  58  56  0.611
Group1  58  58  1
Group1  58  61  0.60113
Group1  58  63  0.534858
Group1  61  56  0.448
Group1  61  58  0.60113
Group1  61  61  1
Group1  61  63  0.622206
Group1  63  56  0.3587
Group1  63  58  0.534858
Group1  63  61  0.622206
Group1  63  63  1
Group2  78  78  1
Group2  78  85  0.622703
Group2  78  94  0.622244
Group2  78  23  0.687826
Group2  85  78  0.622703
Group2  85  85  1
Group2  85  94  0.697313
Group2  85  23  0.612089
Group2  94  78  0.622244
Group2  94  85  0.697313
Group2  94  94  1
Group2  94  23  0.603834
Group2  23  78  0.687826
Group2  23  85  0.612089
Group2  23  94  0.603834
Group2  23  23  1
Group3  27  27  1
Group3  27  80  0.118955
Group3  27  147 0.32038
Group3  27  9   0.335264
Group3  80  27  0.118955
Group3  80  80  1
Group3  80  147 0.430287
Group3  80  9   0.406426
Group3  147 27  0.32038
Group3  147 80  0.430287
Group3  147 147 1
Group3  147 9   0.546452
Group3  9   27  0.335264
Group3  9   80  0.406426
Group3  9   147 0.546452
Group3  9   9   1

My code
corr_matrix = []

for x, group in df.groupby('diabetes_group'):
    matrix = group[["diabetes_group", 'DM-1', 'DM-2', 'corr']]
    matrix_pivot = pd.pivot_table(matrix, values = 'corr', index=['DM-1'], columns=['DM-2'])
    d = matrix_pivot.style.background_gradient(cmap='Set3')

Thanks a lot  M__

corr_matrix = []

for x, group in df.groupby('diabetes_group'):
    matrix = group[["diabetes_group", 'DM-1', 'DM-2', 'corr']]
    matrix_pivot = pd.pivot_table(matrix, values = 'corr', index=['DM-1'], columns=['DM-2'])
    d = matrix_pivot.style.background_gradient(cmap='Set3').set_caption(x)
    display(d)

Sorry I am really late answering, the set_caption was cool.
Background What we are talking about here is the style Class in pandas of which set_background is one method. Usually, it is for outputting a pandas dataframe as a webpage. If I was front-end developer it would be cool. In this example it is used to colour the different bits of the dataframe, presumably to make it easier to read a lot of data.
The thing I would have done is dump to Excel. This will preserve the colouration. The bits contained within # are my codes.
from Pathlib import Path
import pandas as pd

##########
path = Path('/path/todir/')
##########
for x, group in df.groupby('diabetes_group'):
    matrix = group[["diabetes_group", 'DM-1', 'DM-2', 'corr']]
    matrix_pivot = pd.pivot_table(matrix, values = 'corr', index=['DM-1'], columns=['DM-2'])
    d = matrix_pivot.style.background_gradient(cmap='Set3').set_caption(x)

    ##########
    pathout = Path(path, x + '.xlsx')
    d.to_excel(pathout)
    ##########

What was happening originally was the 'loop' was overriding the diabetes Group in previous loop, hence appeared to be disappearing. What will happen now is an excel file will be produced separately for each group and the file name labeled as Group1.xlsx, Group2.xlsx etc ...
Here I dumped my test output to HTML simply via to_html and then to pdf.

Notes
The variable corr_matrix = [] does not take part in this code. It may be important to other bits of code. I'm not sure matrix variable does much, but I've not checked.
The issue with IPython display is the output is not saved, so the code needs to be run each time the analyst wants to examine the correlation coefficients between DM-1 and DM-2.

My personal dataframe was
df = pd.DataFrame({'diabetes_group':['Group1','Group1','Group2','Group2', 'Group3','Group3'],'DM-1':[56,56,56,56,56,56], 'DM-2':[56,58,61,63,56,58], 'corr':[0.2,0.3,0.5,0.3,0.2,0.1]})


text:
How to z-transform Fst and -log(p) values for genome wide selection scan?
I saw in the literature a lot genome wide scans showing $F_{ST}$ values in a 'z-transformed' format which seems to be bounded between 0 to infinity.
Examples are found easily searching for 'genome-wide scan $F_{ST}$' and here is a paper where I could find such z-transformed $F_{ST}$:

Identification of genomic variants putatively targeted by selection during dog domestication (see figure 1)
Genomic divergence during feralization reveals both conserved and distinct mechanisms of parallel weediness evolution (see figure 3)

How are the z-transformed $F_{ST}$ values calculated?
The only reference I could find so far is this: Calculate $F_{ST}$to draw Manhattan plot. And this Population analysis of the Korean native duck using whole-genome sequencing data where they just take the ($F_{ST}$ - mean $F_{ST}$)/sd_$F_{ST}$. But this is giving a lot of negative values, so should this be also surrounded by an absolute value?
Then is it possible to use the z-scores to get -log(p-values)?

The standard score can be estimated by using the sample mean and sample standard deviation as estimates of the population values. Just as you wrote it in the question.
The z-score is given by
$ z = \frac{x - \bar{x}}{s}$
where:
$x$ is the sample observation
$\bar{x}$ is the mean of the sample,
$s$ is the standard deviation of the sample.
Negative values mean that particular observation is $z$ standard deviations below the mean value, and negative $z$-scores are very much expected. The process of standardization is a simple linear transformation and thus does not remove or change any degree of skewness in the original data. I guess your data is skew from the begging to values below the mean.
Yes you can compute the exact $p$-values from the $z$-scores and then apply a negative log transform. You can use a $z$ pre-computed table for this. Check this video for an example.

Thanks for this @JRodrigoF this guy (YouTube link) for teaching critical probability is a teaching maestro.

@JRodrigoF is the correct answer. I'll attempt to explain why the z-score is being used, because it is not immediately obvious.
$F_{ST}$ is not trivial for anyone to understand per se. It is extremely useful for assessing migration - but what is the actual index that Ron Fisher devised?
Simulations suggest that $F_{ST}$ = 0.2 is threshold above which which strong population structure is considered. Thus $F_{ST}$ = 0 panmixia and $F_{ST}$ = 1 is zero migration. How to interpret everything in-between is difficult, except for the 0.2 threshold if the theoretical assumptions in the threshold calculation apply to real-world data. Thus the Z-score is an interesting concept.
In other words, the 0.2 threshold simulation suggests $F_{ST}$ is a non-linear index, which is very interesting in understanding it. Non-linear indices are a bit hairy because we all think they are linear, so 0.1 to 0.3 is one increment of migration, 0.3 to 0.5 another increment of migration. This is not true according to the simulation studies.
Z-scores get around this difficult by simply standardising the combined $F_{ST}$ and thats why they are using it.
It becomes particularly interesting if $F_{ST}$ is calculated per allele because then "fitness traits" can then be tracked. The allele doesn't necessarily need to be the actual fitness trait because of molecular "hitchhiking" (old John Maynard Smith terminology). Thus $F_{ST}$ Z-score within a population within a genome is the beginning of a datamining process. So example, the domestication of dogs why it becomes such a big deal, particularly if the comparison is a dog breed that is closest to the wolf genome. The actually traits selected for in dogs can now be identified. It would also operate between genomes, to identify outliers in dog dispersal, but the behavioural (behavioral) modifications that make dogs dogs is very interesting.

The question is are super small probabilities generated to make log meaningful? I'm not sure in this instance probabilities well beyond very highly, highly significant are generated. If they are then a log transformation is useful.
For background the log-scale is used for transforming very small probabilities rather than use e nomenclature. It makes them more understandable for comparative purposes and allows you to do cool things like add and subtract them instead of multiplying and dividing. The smaller the probability value the higher the transformed value ... so in this case the higher the value the more interesting the result.
However, if the field simply does not use log transformation then it's not cool, because no-one will be able to comparatively interpret the transformed values. In phylogenetics it's done all the time, because the probabilities are really small, but anyone can immediately understand the values - because everyone uses this system. A key part of stats is clarity, not necessarily for those outside the field, but certainly for those within the field.

*very minor quibble it was Sewall Wright not Fisher who invented f-statistics

@user438383 okay, but F is for Fisher (the statistic is named after him), so he must have done something with the development of the stat to get his name recognised like this.

Actually these are not Fisher's F statistics (in fact Fisher didn't develop the F-test, it was just named after him), but Sewall Wright's family of F-statistics, or fixation indices. Confusing, I know...

@user438383, ok I dunno. The reality is Fisher invented so much of classical statistics prior Bayes and maximum likelihood, discovering something he didn't invent really surprises me.


text:
How to rename a file based on a directory name?
I have a lot of FASTA files that I recently scaffolded with RagTag. The output of RagTag includes a directory for each sample and the corresponding FASTA file written as "ragtag.scaffold.fasta"
I would like to go through each file and rename it like so:
Original Path with File name:
directoryName_ID/ragtag.scaffold.fasta

Changed File name:
directoryName_ID/ID_ragtag.scaffold.fasta

As you can see, the directory name includes the ID number I would like to use to rename the fasta file with (along with an underscore to space it out). Is there a way to loop through all of the directories, using that last number after the "_" in the directory name and to use that to rename my fasta file?
Here is what I have tried so far:
for filename in ./directory*/ragtag.scaffold.fasta; do
    mv filename ${./directory%*_*}_ragtag.scaffold.fasta;
done

It doesn't work at all but I know I need to use mv to rename the filename. However, when it comes to referencing the directory name the file is located in I get lost. I also have a .txt file that contains a list of IDs if that might be easier to reference. I'm still kinda new to the command line so if anyone could help that would be greatly appreciated!

What operating system are you running on?

Just execute the following immediately above the directories of interest, (post edited to include OPs "constant term")
for d in *; do
  mv "$d"/ragtag.scaffold.fasta
 "$d/${d/000mergedCluster_/}_ragtag.scaffold.fasta";
done

The key thing I assumed is that the directory name comprises a "constant term" and is suffixed by a "variable term".
Just to note the $d here is directory not file name - because that is constant in each directory.

Comments...
$d/${d/000mergedCluster_/}

The key part of the code is above.
What its doing is saying from the variable $d delete the term 000mergedCluster_ (or whatever the "constant term" is)

Thank you this works great! Only thing is that it includes the "constant term" in the new renamed file...is there a way to change it so only the "variable term" (whats after the _ in the original file name) is added to the file name?

@rimo What is your "constant term" in the directory name?

The whole name is something like "000mergedCluster_1234" so the constant term is "000mergedCluster_" and the variable term is "1234"

@rimo Can you try it now? If you copy and pasted the edited code that should work

She's perfect! Sorry that was silly of me but I completely understand what this command is doing now. Thank you so much for your help!! This works great!

The asterisk wildcard by itself will match all files and directories (except dot files/folders). There's no guarantee that $d is actually a directory. Better to use a more specific pattern, for example: for d in 000mergedCluster_*.

@Steve agreed. That should have been updated, but at the time the directory name was not known, however for d in *_* would have been better.

Assuming the output is only a list of directories and each directory includes a file called "ragtag.scaffold.fasta", you could use brace expansion to generate the SOURCE and DEST strings:
for i in directoryName_* ; do
    echo mv "${i}"/{,"${i#directoryName_}"_}ragtag.scaffold.fasta
done

But the above will only work using Bash and Zsh, and maybe some other shells. Without Bash/Zsh, just use:
fasta="ragtag.scaffold.fasta"
for i in directoryName_* ; do
    echo mv "${i}/${fasta}" "${i}/${i#directoryName_}_${fasta}"
done

Remove the echo to perform the rename.

If you want to use Perl:
#!/usr/bin/perl
use 5.030;
use warnings;

use File::Find::Rule qw( );

for my $file (File::Find::Rule->name("ragtag.scaffold.fasta")->in(".")){
  my $newfile = $file;
  $newfile =~ s{(.+_(\d+)/)(.+.fasta)}{$1$2_$3};

  if (-e $newfile) {
    warn "can't rename $file to $newfile: $newfile exists\n";
  }elsif (rename $file, $newfile) {
    # do nothing
  }else{
    warn "rename $file to $newfile failed: $!\n";
  }
}

Perl one-liner:
ls */filename | perl -lane 'print "mv $1/$3 $1/$2_$3" if /(\S+_([^\/_\s]+))\/(\S+)/'


text:
Is it possible to do homology inference across species using different kinds of NGS data?
Background: I have a list of species that I want to put through homology inference. The goal of homology inference is to investigate the evolution of a trait on a species tree. I want to use the quantitative aspect of the data like CNV and recreate nodal CNV using cafe or count. And use gene trees of the homology to see if I can discern divergence patterns across gene families of interest.
The issue is some species only have RNA-seq generated transcriptions (peptides and cDNA) while others have genomes. The issue of missing introns aside, how can I circumvent this gap?
Should I:

Just use the cDNA data with genomic data? compare nucleic acid to nucleic acid. I need to eventually make gene trees as well and as such MSA is important.The list of species is huge and I anticipate seg-faults during MSA if not during homology inference itself.

Use emboss or something else to convert the nucleic acid data to peptides.And then compare peptides to peptides? Which makes the homology inference → MSA (hopefully pal2nal and gene tree inference later) a little easier. I have never done translation before so I do not know the pitfalls per se.

My list of species is green plants.
Can anyone share their experience if you have done something like this before?

I think you need to use proteins/amino acids.

Personally, I would use blastx against a local protein database. This converts query cDNA->amino acids and therefore combines points 1 and 2 without need to perform a translation in advance.
Genetic distance and saturation might become an issue here depending whether the search is outside the family and especially Order

MAFFT vs Muscle Super5. I would try "Muscle5" or rather the latest Muscle under the "Super5" algorithm. It outperforms everything.

Thanks for suggesting Muscle5, I will have to try it out this time.

this is the right workflow. deleting my answer.

This, or tblastx to compare translated peptide which avoids the need for a local protein database.


text:
pyScenic CLI for ctx is giving error: Not a single module loaded
Hi I ran pyScenic's ctx via the command in command prompt:
pyscenic ctx /home/vboxuser/pyScenic/adj.csv \
    /home/vboxuser/pyScenic/mm10__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather \
    --annotations_fname /home/vboxuser/pyScenic/motifs-v9-nr.mgi-m0.001-o0.0.tbl \
    --expression_mtx_fname /home/vboxuser/pyScenic/csk_filtered_scenic.loom \
    --mode dask_multiprocessing \
    --output reg.csv \
    --mask_dropouts \
    --num_workers 20

It gives me a bunch of warnings and then an error, as shown below:
2023-01-17 09:54:52,135 - pyscenic.cli.pyscenic - INFO - Creating modules.

2023-01-17 09:54:52,165 - pyscenic.cli.pyscenic - INFO - Loading expression matrix.

2023-01-17 09:54:55,458 - pyscenic.utils - INFO - Calculating Pearson correlations.

2023-01-17 09:54:55,482 - pyscenic.utils - WARNING - Note on correlation calculation: the default behaviour for calculating the correlations has changed after pySCENIC verion 0.9.16. Previously, the default was to calculate the correlation between a TF and target gene using only cells with non-zero expression values (mask_dropouts=True). The current default is now to use all cells to match the behavior of the R verision of SCENIC. The original settings can be retained by setting 'rho_mask_dropouts=True' in the modules_from_adjacencies function, or '--mask_dropouts' from the CLI.
    Dropout masking is currently set to [True].

2023-01-17 09:54:56,273 - pyscenic.utils - INFO - Creating modules.
/home/vboxuser/mambaforge/lib/python3.10/site-packages/pyscenic/utils.py:244: FutureWarning: Not prepending group keys to the result index of transform-like apply. In the future, the group keys will be included in the index, regardless of whether the applied function returns a like-indexed object.
To preserve the previous behavior, use

    >>> .groupby(..., group_keys=False)

To adopt the future behavior and silence this warning, use

    >>> .groupby(..., group_keys=True)
  df = adjacencies.groupby(by=COLUMN_NAME_TARGET).apply(
/home/vboxuser/mambaforge/lib/python3.10/site-packages/pyscenic/utils.py:244: FutureWarning: Not prepending group keys to the result index of transform-like apply. In the future, the group keys will be included in the index, regardless of whether the applied function returns a like-indexed object.
To preserve the previous behavior, use

    >>> .groupby(..., group_keys=False)

To adopt the future behavior and silence this warning, use

    >>> .groupby(..., group_keys=True)
  df = adjacencies.groupby(by=COLUMN_NAME_TARGET).apply(
/home/vboxuser/mambaforge/lib/python3.10/site-packages/pyscenic/utils.py:244: FutureWarning: Not prepending group keys to the result index of transform-like apply. In the future, the group keys will be included in the index, regardless of whether the applied function returns a like-indexed object.
To preserve the previous behavior, use

    >>> .groupby(..., group_keys=False)

To adopt the future behavior and silence this warning, use

    >>> .groupby(..., group_keys=True)
  df = adjacencies.groupby(by=COLUMN_NAME_TARGET).apply(

2023-01-17 09:54:56,350 - pyscenic.cli.pyscenic - ERROR - Not a single module loaded

I tried uninstalling and reinstalling pyScenic, but it didn't change anything. It shouldn't be pyScenic overall either that is wrong, since I ran pyScenic grnboost command and that worked fine.
It seems like it might be specific to the pyScenic commands. Has anyone run into the same issues?
Thanks!

Background information is here https://pyscenic.readthedocs.io/en/latest/tutorial.html


text:
How can I predict the folding of a synthetic protein?
I want to predict the folding of a chimeric antigen receptor sequence. These sequences are completely made out of small parts of different proteins. Do you know if AlphaFold can be used for this?

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.

From @MaximillianPress's response, the answer is no because a chimera has no history the domains evolving together and that is a central part of the algorithms assumptions. Please do remember to mark an answer as "accepted" when the question has been successfully answered.

I'd consult references about AlphaFold. For example, wikipedia has a "Limitations" page on the tool, though the page doesn't mention anything about synthetic proteins.
The hesitation that I would have is that the tool relies to some extent on coevolutionary models of amino acids to infer interactions at both training and prediction steps, so throwing a synthetic protein at it may confuse it. See e.g. slide 10 here.
In an interview, they acknowledge this rather directly:

Now, if you have a protein which is completely new, which does not have an evolutionary history—a synthetic protein, which has just been designed— then our system would not be as accurate as it would be for existing proteins.

I am not totally clear on whether this is due to learning a latent representation of constraints on similar proteins, or whether the tool attempts to directly find homologs. The paper methods suggest to me that it is the former, but I'm not sure. They describe training on only chunks of a big MSA. If that's the case, I'd guess that they still can use information about the natural chunks of your synthetic protein, and aren't expecting to be able to do a full-length protein alignment.
So I think that you can expect the prediction to be relatively low quality. It is still probably much better than any alternative way of predicting structure!
As always, not a bad idea to post an issue on github asking this question.

This is a very good point indeed, it will rely on convolution as part of its algorithm and chimeras will not follow this assumption.


text:
How to fix the Error in Kruskal Wallis rank sum test for comparing the guides in python? Also to check whether using the right test for my data?
df
gene_name   guide_1 guide_1_new Correlation
MMP-1   A   A   1
MMP-1   A   B   0.426115
MMP-1   A   C   0.522499
MMP-1   A   D   0.431587
MMP-1   B   A   0.426115
MMP-1   B   B   1
MMP-1   B   C   0.60113
MMP-1   B   D   0.534858
MMP-1   C   A   0.522499
MMP-1   C   B   0.60113
MMP-1   C   C   1
MMP-1   C   D   0.622206
MMP-1   D   A   0.431587
MMP-1   D   B   0.534858
MMP-1   D   C   0.622206
MMP-1   D   D   1

I used four guides in an experiment (data = df). I calculated correlation between them based on the their log fold changes. I want to check which guide is outperforming or underperforming in comparison to other three guides in my experiment. I tried Kruskal and Wallis test to do this analysis with the script given below. However, I got an error ValueError: Samples must be one-dimensional. Also is this the right test to do this analysis?
My script
from scipy import stats

x = df[['guide_1']]
y = df[['guide_1_new']]
z = df[['Correlation']]
stats.kruskal(x, y, z)

I am a little confused, you are passing the columns of guide_1 and guide_1_new as samples to scipy.stats.kruskal? Those are non-numeric data. And the correlations, while numeric, suggest to me that you are doing something a little different from what Kruskal-Wallis is implementing. It looks like you might be trying to pass in the categories as vectors of strings, but I don't think it works that way (?). And I don't think applying it to correlation coefficients will be very interpretable. Not an expert on scipy.stats but this doesn't look right. I'd suggest going back to raw data.

Answered below, if you are stuck I can talk you through the calculation if you perform a df.to_dict() and post the output.

Reformatting the data for the KW and doing the test is easy.
Using "guide_1" alone and pivot() to get A = [1, 0.426115 ...
df2 = df.pivot(columns='guide_1', values='Correlation')

This will now present A, B, C, D etc ... into columns...
x = df2[['A']].dropna()
y = df2[['B']].dropna()
z = df2[['C']].dropna()
a = df2[['D']].dropna()

stats.kruskal(x, y, z, a)

That will work ... whether that is a good test is a complicated question because this is log transformed data. Transformation is used in parametric tests. Removing the log transformation is strongly recommended because KW works with raw data.
You can switch the calculation to guide_1_new instead. It will give a different answer ... but I don't know the experiment. The original guide looks cool.
Looking at the data - the test should work because I believe KW is 'unpaired' (you need to check) and looking at the value the distribution looks pretty similar between A, B, C and D. However, just because the KW says 'yes' doesn't mean thats a robust answer (log data). If the untransformed data says 'yes' ... thats a robust answer for KW. Just check the 'unpaired' stuff - thats important for the result.

Personal note, I originally used a third df but was not needed,
df3 = pd.concat([df2[col].sort_values().reset_index(drop=True) for col in df2], axis=1, ignore_index=True)

I'm a little confused here, it looks to me in the OP like the numeric input is correlation coefficients from comparing A to A, A to B, etc. So for A to A, you get an unsurprising correlation of 1, which is now included in the A group, which doesn't seem like a very informative input. Does OP really care about whether guide_1 group A is more correlated to other groups than those groups are among themselves (which is what this is testing for, effectively)? Doesn't make sense to me, shouldn't OP go back to raw (pre-correlation) data? Probably missing something here.

@MaximilianPress KW will give a result representing the data. Agreed non-parametrics prefer raw data. The presence of A->A , B->B will of course bias the test, but the gist is A ~ B ~ C ~ D - and will not be 'thrown out of court'. Capturing log-change or non-log change isn't easy. A 2-way ANOVA type (4 columns vs 2 rows) would do this - absolutely, but its complicated and non-parametrics are not good at this calculation (Friedman test, statisticians avoid it). Using 1-way ANOVA on this specific data would be a mistake. Getting it 'spot on' is an investment. Right now 'it's ok'.


text:
How to highlight the specific peptide sequences after performing multiple sequence alignment for the fasta file?
I have the peptide sequences and fasta files separately. I first aligned the fasta files using msa package. After that I'm trying to highlight the peptide sequences in the multiple sequence alignment output. I couldn't find a way to do that. Any suggestions on how to do this would be very helpful.
My R script so far:
# load required packages
library(seqinr)
library(msa)
library(dplyr)

# read the file containing the list of peptide sequences
peptides <- read.table("./peptides/AB_hydrolase-1_domain-containing_protein.txt", header=F, stringsAsFactors=F)

# read the file containing the list of fasta files
#hemoSeq <- readAAStringSet(system.file("examples/HemoglobinAA.fasta",package="msa"))
fasta_file <- "./fasta/AB_hydrolase-1_domain-containing_protein.txt"
fasta_data <- readAAStringSet(fasta_file)

# multiple sequence alignment with ClustalW algorithm
msa <- msa(fasta_data, method = "ClustalW")
```

Hi and welcome. Could you please define what you "specific peptide sequences"? Perhaps also include the title of the question in the question.

@M__ I have some short peptide sequences belonging to these fasta sequences, i want to highlight these peptides positions to know their locations in the msa output

Checkout an MSA visualisation and analysis package that you can load the alignment into. It may not be the newest tool, but I've used Jalview (https://www.jalview.org/) in the past

I have used jalview myself. But I want to do it in a R program


text:
How to remove empty files from a directory?
I have a directory of FASTA files, some of which are empty (placeholders that were output from a pipeline) and I would like to delete the empty files using a for loop on the command line.
From what I have seen so far I can use the wc -l [file] approach to see "if there are lines present, delete the file" but I want to make sure the syntax is right before deleting important files. I was thinking it might be easier for the command to look for the " > " symbol found at the header in FASTA files and if it's not present delete the file.
Here is what I have so far:
for file in *.fasta; do "if wc -l < 1 rm file"; done

Or
for file in *.fasta; do "grep > | rm file"; done

You'll want to use the find command, as the answer from @Ram RS shows. But on a separate note, neither of the above commands are syntactically correct. For a safe way to use these commands & test them, substitute ls for rm, to simply list files rather than delete them.

If you try find with -links and let me know

What operating system are you using? Not all find implementations have the -empty or -delete flags.

This isn't the most efficient approach, using find will be easier and faster as long as your find supports the -empty option, but since this is what you asked for, here is how to do it based on the number of lines:
for file in *.fasta; do
 [[ $(wc -l < "$file") -eq 0 ]] && echo rm -- "$file"
done

If you want to recursively descend into subdirectories, use:
shopt -s globstar
for file in  *.fasta */**.fasta; do.
  [[ $(wc -l < "$file") -eq 0 ]] && echo rm -- "$file"
done

If that prints the right files, remove the echo to actually delete them.
Note that this will also delete files that are not empty but don't have a newline character. These are not valid "text files", but strictly speaking they aren't empty, they just have no lines as lines are defined by the \n character. For example, this file would be deleted:
printf 'aaa' > file

The potential problem with this approach in Nextflow is there are "extensive" (or more than I think reasonable) numbers of subdirectories in the output and the names get complicated, which precisely is why I avoid it  (I'd rather construct my own Dockerfiles). find will run straight through the directory structure and the depth can be set. Thats why in my opinion it's the favoured solution here.

@M__ yes, absolutely. If the OP's find supports the -empty and -delete, it is the best tool for the job. You could easily make this recursive with shopt -s globstar; for file in  *.fasta */**.fasta; do....

thank you so much for this line of code!! and the comments really help too. I am running this with a mac OS which is why I think the find function doesn't work for me. It definitely seems like the easier option so I might try and find a workaround with an aws instance in the future but for now this works like a charm. Thank you!

@rimo yes, exactly, macOS ships with a different find implementation that might not support all of the flags used like -maxdepth. -empty and -delete.

Always Google before asking on forums. I googled "delete empty files linux" and this was the second result: https://www.baeldung.com/linux/delete-empty-files-dirs
# Print and delete empty files just in the current directory
$ find . -maxdepth 1 -type f -empty -print -delete
# Print and delete empty files in the current directory and all sub-directories
$ find . -type f -empty -print -delete

I saw that link too and tried it but it didn't work for some reason

What happened when you tried it?

it doesn't do anything...nothing output in the STDOUT and nothing changed in the directory

Are you running the command from the right directory? If it does not delete anything, but you're sure you have "empty" files, check if the files have blank spaces or some such invisible characters. cat -A them.

@rimo Ram's solution only works on Linux systems, what are you using?

I'm using a mac OS which is probably why that doesn't work for me

It should work on macOS too - it's not because of the OS, it's because of the file type. Empty aliases are not empty files, much like how empty directories are not empty files. Define your tasks more specifically and you'll be able to get to an accurate answer quicker.

As already stated, one way is to use the GNU find (the default on Linux) command. This will find and print all empty files in the current directory:
find . -maxdepth 1 -type f -empty -print

Append the -delete action to the command above when you're ready to remove them. If you need to find (and then delete) empty files recursively, simply remove the -maxdepth 1 option.

In a comment, you mention that these are actually "alias" files which are the results of a Nextflow pipeline. Unless the default publishDir mode has been changed, Nextflow will publish output files to the publishDir using symbolic links to the files in the working directory (i.e. ./work). I suspect what you are trying to do is to find and remove the links to these empty files, but not the files themselves. In this case, I think what you want is:
find . -type l -exec test ! -s '{}' \; -print

Again, append the -delete action to the command above to actually remove the symlinks. If you instead need to delete files from Nextflow's working directory, use nextflow clean in conjunction with nextflow log to safely clean up the failed runs.

@Ram RS has identified the solution here via find, if that does not work then the following will work,
find . -maxdepth 1 -type f -empty | xargs rm

Edit 2
The solution is obvious ...
find . -type l -delete

The code under Edit is is better if the OP knows in advance precisely how many links the symbolic link file has, basically its always one. However, we don't know, thus -type l will delete every symbolic link.

To address @Steve's point of inadvertent deletion in case of a null return ... The code was checked prior release and this didn't occur (and also now rechecked).

I think the issue I'm facing is that these files are "Alias" files which are the results from a Nextflow script I wrote. Therefore, the file has size (~180 bytes) but is ultimately empty inside. I think the solution I am looking for is more specific than using the "find" command because it's not working for me. I would like to try and use a solution that actively looks in the files and sees whether there is text there (i.e. wc -l < 1) and if not, delete the file.

Or change the scope of what you are doing to only make softlinks to files that are greater than 0 bytes in size, rather than deleting the ones you already have made.

If find doesn't return any results, a pipe to xargs rm -rf could remove everything in your current directory.

@Steve the point was addressed above

if the -delete fails, that is most likely because the OP is using a find version that doesn't have it, so will very likely also not support the -empty. Also, note that this will fail on weird file names (which shouldn't be an issue with fasta files, but you never know). More importantly, your "edit" section does not remove symbolic links, the -links predicate is about the number of hard links. Since most inodes only have a single hard link (the file), your command will simply delete all files in the directory.

@M__ no, it will not remove symbolic links. Or, more accurately, it will remove all files, which includes symbolic links so in that sense it will remove symbolic links. See man find where the -links option is documented and where it explains that it refers to hard links and also just try touch file{1..10} and then run find . -type f -links 1 -print -delete and you will see that all 10 of the files you created will be removed. And yes, the first point applies to all, but the OP commented here hence I left the comment here to explain why it failed for them.

Okay I err ....

I don't understand how an Alias file (or even a symlink for that matter) can be considered empty. Is the alias/link empty when the file pointed to is empty? Or is it empty when the file pointed to does not exist? I'm unable to understand OP's situation.


text:
How to write FASTA records using "Bio.SeqIO.write()"
Update:
Biopython document says that "Bio.SeqIO.FastaIO.FastaWriter" class is obsolete. Now my question becomes how to I use Bio.SeqIO.write() fucntion to achive the same thing I asked.
I use the following code to extract specific FASTA records from a file based on a list of FASTA IDs given in another file. The script works just fine. However, I want to able to write the results using "Bio.SeqIO.FastaIO.FastaWriter" class into a file instead of Python's print or write function. How can I do that in my last if statement?
from Bio import SeqIO
#import pprint
import textwrap

my_id_file     = open('sample_files2/fasta_id_records.txt','r')
my_fasta_file  = open('sample_files2/fasta_file.fasta','r')

my_dictionary = {} # fasta IDs are keys, value can be anything.
for line in my_id_file:
    my_dictionary[line[:-1]] = 'value'
#pprint.pprint(my_dictionary)

for seq_record in SeqIO.parse(my_fasta_file, "fasta"):
    if seq_record.id in my_dictionary:
        sequence = str(seq_record.seq)
        fasta_record = textwrap.fill(sequence, width=60)
        print(f">{seq_record.id}\n{fasta_record}")

my_fasta_file.close()
my_id_file.close()

I was considering modifying the title. I am glad you edited it @gringer. Thank you.

I figured it out:
from Bio import SeqIO
#import pprint
import textwrap

my_id_file     = open('sample_files2/fasta_id_records.txt','r')
my_fasta_file  = open('sample_files2/fasta_file.fasta','r')
result_file    = open('sample_files2/results.fasta', "w")

my_dictionary = {} # fasta IDs are keys, value can be anything.
for line in my_id_file:
    my_dictionary[line[:-1]] = 'value'

#pprint.pprint(my_dictionary)
for seq_record in SeqIO.parse(my_fasta_file, "fasta"):
    if seq_record.id in my_dictionary:
        #use this line to write standard fasta record using SeqIO.write
        SeqIO.write(seq_record, result_file, "fasta")

my_fasta_file.close()
my_id_file.close()
result_file.close()

I'm not 100% sure you can because that will be taken care within SeqIO.write(). The obvious seq_record.seq = textwrap.fill(str(seq_record.seq) , width=60) is unlikely to work.

The only way to do it is within a straight write to disk in you first code ... but (obviously) write rather than print.

it wraps to 60 by default in Biopython 1.8 at least, if saves as 'fasta-2line' uses a single line for the seq

You could instead pass in a generator that yields SeqRecords. This avoids storing the entire list of  sequences in memory and is usually the way to go assuming that you need to handle large files containing thousands of records. This example use a generator expression, which is a bit like a list comprehension, but instead of creating a list, it creates a generator:
import argparse
import contextlib
import pathlib
import sys

from Bio import SeqIO

def get_argument_parser():

    parser = argparse.ArgumentParser()

    parser.add_argument(
        "fasta_file",
        type=pathlib.Path,
        metavar='FILE',
        help="Extract sequence records from FILE",
    )
    parser.add_argument(
        "ids_file",
        type=pathlib.Path,
        metavar='FILE',
        help="Sequence record IDs/names to extract, one per line",
    )
    parser.add_argument(
        "-o",
        "--output",
        type=pathlib.Path,
        metavar='FILE',
        default="-",
        help="Write extracted sequences to FILE (default: stdout)",
    )
    parser.add_argument(
        "-f",
        "--force",
        action='store_true',
        help="Overwrite the output file if it exists",
    )

    return parser

def main():

    parser = get_argument_parser()
    args = parser.parse_args()

    with contextlib.ExitStack() as stack:

        ids_file = stack.enter_context(open(args.ids_file))
        fasta_file = stack.enter_context(open(args.fasta_file))

        if args.output != '-':
            out = stack.enter_context(open(args.output, 'w' if args.force else 'x'))
        else:
            out = sys.stdout

        record_ids = {line.rstrip('\n') for line in ids_file}

        fasta = SeqIO.parse(fasta_file, "fasta")
        records = (rec for rec in fasta if rec.id in record_ids)

        SeqIO.write(records, out, "fasta")

if __name__ == "__main__":
    main()

cool never eard about ExitStack till now https://docs.python.org/3/library/contextlib.html#contextlib.ExitStack

A better description about what it does is here https://www.rath.org/on-the-beauty-of-pythons-exitstack.html @pippo1980

I think I need to ask another question on this code. It's certainly professional code. What I'm wondering is how the if-condition tuple comprehension is a generator because thats the "workhorse" of the code.

@M__ The parentheses just make it a generator expression. Brackets would have made it a list comprehension: https://peps.python.org/pep-0289/

@Steve, gave me a code to study.

my attempt at changing the way Biopython formats fasta file sequence width,
read into code to figure out standard behaviour:
from Bio import SeqIO, __version__
#import pprint
# import textwrap

print('Biopython version : ' , __version__)

# my_id_file     = open('sample_files2/fasta_id_records.txt','r')
# my_fasta_file  = open('sample_files2/fasta_file.fasta','r')
#result_file    = open('sample_files2/results.fasta', "w")
### CHANGED TO:

my_id_file     = open('fasta_id_records.txt','r')
my_fasta_file  = open('fasta_file.fasta','r')
result_file    = open('results.fasta', "w")

my_dictionary = {} # fasta IDs are keys, value can be anything.
for line in my_id_file:
    my_dictionary[line[:-1]] = 'value'

#expected lenght of sequences in my fasta
length = 10

#pprint.pprint(my_dictionary)
for seq_record in SeqIO.parse(my_fasta_file, "fasta"):
    if seq_record.id in my_dictionary:

        ########################################################################################
        #use this line to write standard fasta record using SeqIO.write

        # SeqIO.write(seq_record, result_file, "fasta-2line") # writes fasta in 2 line format

        # SeqIO.write(seq_record, result_file, "fasta") # writes fasta with width of 60 as standard biopython format

        ### NOT POSSIBLE !!!!!!!!!!!!
        # seq_record.seq = textwrap.fill(seq_record.seq , width=10) ### --> AttributeError: 'Seq' object has no attribute 'expandtabs'

        # print(seq_record.seq, type(seq_record.seq), dir(seq_record.seq))

        sequence = str(seq_record.seq)

        result_file.write('>'+seq_record.id+'\n')

        while len(sequence) > 0:
            result_file.write(sequence[:length]+'\n')
            sequence = sequence[length:]
        result_file.write('\n')

        ### if you try to add EOL inside seq_record.seq of course you get
        ### TypeError: SeqRecord (id=seq1) has an invalid sequence.
        # SeqIO.write(seq_record, result_file, "fasta") # writes fasta with sequence width as 60   !!!! TypeError: SeqRecord (id=seq1) has an invalid sequence.

my_fasta_file.close()
my_id_file.close()
result_file.close()

input files :
fasta_id_records.txt:
seq1
seq2
seq3

 fasta_file.fasta:
>seq0
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
>seq1
TTTT
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
TTTTTTTTTTTTTTT
>seq3
CCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCGGGGGGGG
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG

output file:
v :
>seq1
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTT

>seq3
CCCCCCCCCC
CCCCCCCGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGCCCCC
CCCCCCCCCC
CCCCCCCCCC
CCCCCCCCCC
CCCCCCCCCG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
G


text:
How to get single histogram plot for all groups in different colours?
I am using the below script to get the histogram plot for all values together. However, I want a single plot but different colours for different groups (A, B, C)
data file
Group_name  code-1    code-2    Expression
GroupA  6  6  1.21
GroupA  6  8  0.611
GroupA  6  1  0.448
GroupA  6  3  0.3587
GroupB  8  6  0.611
GroupB  8  8  2.11
GroupB  8  1  0.60113
GroupB  8  3  0.534858
GroupC  1  6  0.448
GroupC  1  8  0.60113
GroupC  1  1  2.2
GroupC  1  3  0.622206

Script
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams["figure.figsize"] = (12,8)

g1 = sns.histplot(df.Expression)
g1.set(
    xlim=(-2,2),
    title="Expression"
)
g1.axvline(0.0, c='black', linestyle="--")
g1.axvline(df.Expression.median(), c='red', linestyle="--")

```

Just change your g1 = sns.histplot(df.Expression)
to g1 = sns.histplot(data=df, x="Expression", hue="Group_name")
Here's the full code:
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd

plt.rcParams["figure.figsize"] = (12,8)

df = pd.DataFrame({
    "Group_name": ["GroupA", "GroupA", "GroupA", "GroupA", "GroupB", "GroupB", "GroupB", "GroupB", "GroupC", "GroupC", "GroupC", "GroupC"],
    "code-1": [6, 6, 6, 6, 8, 8, 8, 8, 1, 1, 1, 1],
    "code-2": [6, 8, 1, 3, 6, 8, 1, 3, 6, 8, 1, 3],
    "Expression": [1.21, 0.611, 0.448, 0.3587, 0.611, 2.11, 0.60113, 0.534858, 0.448, 0.60113, 2.2, 0.622206]
})

g1 = sns.histplot(data=df, x="Expression", hue="Group_name")
g1.set(
    xlim=(-2,2),
    title="Expression"
)
g1.axvline(0.0, c='black', linestyle="--")
g1.axvline(df.Expression.median(), c='red', linestyle="--")

g1.figure.savefig("hist.png")

'''

Thanks, very impressive Seaborn coding


text:
Is there a tool to visualize RNA secondary structure?
I want to manually fold and pair two sequences (U4 and U6). And I want to make niceish figures.
I haven't found a single program that does that.
thanks

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.

Hi @user3392394 and welcome. The question is difficult to understand. Are you trying to perform two individual RNA structures, or is this two sequences to make a single structure? There are lots of programs for RNA structure, what I suspect is you have a specific project in mind, which is not being clearly communicated. Alternatively are you just looking for a RNA 2D structure viewer?

The RNA viewer that is particularly good is VA RNA. This stands for "Visualisation Applet for RNA". Its a java app and run it via,
java -jar /Route2Progam/vaRNA/VARNAv3-93-API.jar

Where Route2Progam is the path to where the jar files is stored, all other info in the run code is literal.
Very nice pictures, point and clicky, drag structures into different shapes and easy to colour ... so for a GUI its great.
They recommend PsuedoViewer for psuedoknots, but if that doesn't mean much ... it doesn't matter.

You could try mxfold2.
https://github.com/mxfold/mxfold2
% mxfold2 predict test.fa
>DS4440
GGAUGGAUGUCUGAGCGGUUGAAAGAGUCGGUCUUGAAAACCGAAGUAUUGAUAGGAAUACCGGGGGUUCGAAUCCCUCUCCAUCCG
(((((((........(((((..((((.....))))...)))))...................(((((.......)))))))))))). (24.8)


text:
Batch correction in differential expression analysis
I have sent two sets (two batches in matter of sending for sequencing) of different samples (plasma) to small RNA-seq to Qiagen company
This is how my meta data look
> head(meta)
           batch group    sex
30678.009 batch2     C   Male
30678.010 batch2     C Female
30941.001 batch1     C Female
30941.002 batch1     C   Male
30941.003 batch1     C   Male
30941.004 batch1     C   Male
> tail(meta)
           batch group    sex
30941.016 batch1    LT Female
30941.017 batch1    LT   Male
30941.019 batch1    LT   Male
30941.022 batch1    LT Female
30941.023 batch1    LT   Male
30941.024 batch1    LT Female
>

I want to extract differentially expressed RNAs between my two groups without batch effect
I am doing this
Is this correct please?
dds <- DESeqDataSetFromMatrix(countData = count,
                              colData = meta,
                              design= ~ batch + group)

dds <- DESeq(dds, test="LRT", reduced=~batch)

res <- results(dds)
resultsNames(dds)

[1] "Intercept"              "batch_batch2_vs_batch1" "group_LT_vs_C"
>

If both batches contain replicates of both groups (LT and C) then the presented code is what you should do, exactly pointed out in the vignette for the LRT: http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#likelihood-ratio-test

Thank you. Please, what do you mean by "If both batches contain replicates of both groups (LT and C)" ? I know each sample is one independent person and we have sent samples to Qiagen in two times

It means exactly like it reads. Both batches need to contain replicates (preferably in comparable numbers) of both groups. Use table(meta$batch, meta$group) to see the distribution contingency table.


text:
How to identify a low proportion cell subpopulation in the single-cell RNA-seq data?
The cell subpopulation that I am interested in only accounts for around 1.2% of the total cells. I have previous FACS experiments that sort out the subpopulation from the samples (using markers CD166+ CD49fhi CD104- Lin-, see images). May I ask if there are any ways to identify this cell population from single-cell RNA-seq data of the tumor epithelial samples using the marker information?

Run the usual steps including clustering and visualization via something like UMAP and then color the UMAP by these four markers. That assumes that the surface marker separation you see in flow holds true on a transcriptional level which is not necessarily the case. For example, hematopoietic surface-defined populations such as the Granulocyte-Macrophage Progenitor (GMP) are well-separated in flow but on transcriptional single-cell cell level span many clusters.
Alternatively, check the tool UCell which can calculate a score per cell based on a combination of positive and negative markers. You can then color the UMAP by this to see whether a single cluster stands out or whether this celltype contributes to many clusters ir is not even present at all. Include some marker combinations that you know are oresent and some nonsense combination that must be absent to have positive and negative controls. That will give you a feel whether that score for your cells is reliable.
This all assumes that the cells are even in the dataset. Depending on cell size and frequency (underrepresentation) they might even have dropped out during library prep.


text:
Exon number between GTF file and ENSEMBL site does not match
Based on the ENSEMBL GTF file, I counted the number of exons per gene as follows:
awk '$3 == "exon" {print}' Homo_sapiens.GRCh38.108.gtf | gffread -F --keep-exon-attrs --table "gene_id","transcript_id",@numexons | sort -k1,1 | datamash -sg1 max 3 > gene_lists/exons_per_gene.txt

An extract of the output is:
ENSG00000034677 11
ENSG00000034693 12
ENSG00000034713 4
ENSG00000034971 4

However, when I check the ENSG00000034713 ENSEMBL site, I count 5 exons. Why is this different?

Thank you all for the clarification. I was looking at the gene sequence section (https://www.ensembl.org/Homo_sapiens/Gene/Sequence?db=core;g=ENSG00000034713;r=16:75566375-75577881) where exons are marked in red (there are 5 red regions for the example gene), but not all transcripts will use all exons as said above. Then the number I have makes sense. Thanks again.

Answer from @terdon, @Ben_Ensembl, converted from comments:
Genes don't have exons, transcripts have exons. There are five transcripts for ENSG00000034713, one with one exon, two with three exons, and two with four exons. No transcripts with five exons though.
The count you retrieved from the GTF will be the number of unique exons for all transcripts of each gene (some unique to a single transcript, some common between two or more transcripts).
[please improve this answer if you can by editing it]


text:
analysing node connectivity difference in a continuous scale (edge weighted by correlations)
I have an edge list object weighted by a correlation value. I would like to know the changes in the node's connectivity on a continuous scale.
In other words, how does the node connectivity differ from the edge weight scale (-1) to the edge weight (+1)? I would like to know how to do this or whether any R package or statistical model could solve the problems etc.
I would like to know whether node connectivity increases or decreases in relation to the correlation weight. I.e. how it looks at the correlation between node degree and edge weight. But I don't know how to do it.
A example edge list object is created in the code below.
data <- structure(list(from = c(5L, 5L, 5L, 1L, 1L, 1L, 1L, 4L, 4L, 4L,
                                2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L),
                       to = c(1L, 2L, 3L, 5L, 4L, 2L, 3L, 1L, 2L, 3L, 5L, 1L, 4L, 3L, 5L, 1L, 4L, 2L),
                   weight = c(runif(18, -1,1))),
                  row.names = c(NA, -18L),
                  class = "data.frame")

I appreciate your help!!
Best,
Amare

if I understand correctly, something like, "do nodes with higher UNWEIGHTED degree also tend to have higher weight on their edges?"


text:
expansions and contractions from OrthoFinder
The input file looks like this, and the complete file can be found here:
OG0000008 NbL07g11380.1 NbL19g07810.1 NbL19g09170.1 NbL19g19070.1 NbQ01g01670.1 NbQ01g03330.1 NbQ01g04070.1 NbQ01g04670.1 NbQ01g05120.1 NbQ01g05870.1 NbQ01g06940.1 NbQ01g07580.1 NbQ01g08860.1 NbQ01g10050.1 NbQ01g10360.1 NbQ01g14200.1 NbQ01g14790.1 NbQ01g16080.1 NbQ01g17760.1 NbQ01g19270.1 NbQ01g19310.1 NbQ01g19390.1 NbQ01g21260.1 NbQ01g21330.1 NbQ01g21740.1 NbQ01g21910.1 NbQ01g23100.1 NbQ01g24620.1 NbQ01g25340.1 NbQ01g26060.1 NbQ01g26320.1 NbQ02g00750.1 NbQ02g03100.1 NbQ02g03420.1 NbQ02g03610.1 NbQ02g03680.1 NbQ02g05120.1 NbQ02g07460.1 NbQ02g08170.1 NbQ02g08330.1 NbQ02g09220.1 NbQ02g09400.1 NbQ02g10620.1 NbQ02g11310.1 NbQ02g14330.1 NbQ02g14460.1 NbQ02g14520.1 NbQ02g15320.1 NbQ02g17090.1 NbQ02g17130.1 NbQ02g20290.1 NbQ02g23070.1 NbQ02g23420.1 NbQ02g24450.1 NbQ02g24480.1 NbQ02g26700.1 NbQ03g00830.1 NbQ03g01970.1 NbQ03g04460.1 NbQ03g06900.1 NbQ03g09530.1 NbQ03g10620.1 NbQ03g12760.1 NbQ03g13450.1 NbQ03g15540.1 NbQ03g15640.1 NbQ03g17180.1 NbQ03g20740.1 NbQ03g21510.1 NbQ03g24670.1 NbQ04g01350.1 NbQ04g01720.1 NbQ04g08420.1 NbQ04g09090.1 NbQ04g10450.1 NbQ04g11470.1 NbQ04g12120.1 NbQ04g14130.1 NbQ04g15440.1 NbQ04g15860.1 NbQ04g16450.1 NbQ04g16620.1 NbQ04g17760.1 NbQ04g19040.1 NbQ04g20020.1 NbQ05g03320.1 NbQ05g04660.1 NbQ05g05970.1 NbQ05g07500.1 NbQ05g08900.1 NbQ05g09760.1 NbQ05g10830.1 NbQ05g11150.1 NbQ05g11340.1 NbQ05g11510.1 NbQ05g11530.1 NbQ05g11780.1 NbQ05g16980.1 NbQ05g18190.1 NbQ05g21710.1 NbQ05g23400.1 NbQ06g01110.1 NbQ06g01430.1 NbQ06g04200.1 NbQ06g04440.1 NbQ06g05330.1 NbQ06g05770.1 NbQ06g05820.1 NbQ06g06700.1 NbQ06g08620.1 NbQ06g09190.1 NbQ06g10460.1 NbQ06g15220.1 NbQ06g15330.1 NbQ06g15700.1 NbQ06g16320.1 NbQ06g16590.1 NbQ06g17590.1 NbQ06g17670.1 NbQ06g20050.1 NbQ07g01030.1 NbQ07g02010.1 NbQ07g04350.1 NbQ07g04900.1 NbQ07g05610.1 NbQ07g06200.1 NbQ07g07110.1 NbQ07g07690.1 NbQ07g08640.1 NbQ07g10390.1 NbQ07g11920.1 NbQ07g14130.1 NbQ07g15590.1 NbQ07g15620.1 NbQ07g16910.1 NbQ07g17130.1 NbQ07g17950.1 NbQ08g00060.1 NbQ08g02240.1 NbQ08g02300.1 NbQ08g02310.1 NbQ08g03290.1 NbQ08g05330.1 NbQ08g09280.1 NbQ08g14890.1 NbQ08g15820.1 NbQ08g15950.1 NbQ08g19830.1 NbQ08g20150.1 NbQ08g22050.1 NbQ08g22620.1 NbQ09g02100.1 NbQ09g02620.1 NbQ09g03950.1 NbQ09g04200.1 NbQ09g06040.1 NbQ09g06640.1 NbQ09g08160.1 NbQ09g08330.1 NbQ09g09660.1 NbQ09g11220.1 NbQ09g13860.1 NbQ09g15180.1 NbQ09g15310.1 NbQ09g16530.1 NbQ09g17900.1 NbQ09g18100.1 NbQ09g18720.1 NbQ09g19280.1 NbQ09g21840.1 NbQ10g00480.1 NbQ10g01350.1 NbQ10g02870.1 NbQ10g03640.1 NbQ10g03730.1 NbQ10g08070.1 NbQ10g09510.1 NbQ10g11010.1 NbQ10g11760.1 NbQ10g12050.1 NbQ10g12060.1 NbQ10g12910.1 NbQ10g19200.1 NbQ10g19930.1 NbQ10g20390.1 NbQ10g20730.1 NbQ10g21080.1 NbQ10g21140.1 NbQ10g24010.1 NbQ11g00310.1 NbQ11g01210.1 NbQ11g01370.1 NbQ11g04610.1 NbQ11g04800.1 NbQ11g06060.1 NbQ11g07820.1 NbQ11g08390.1 NbQ11g09100.1 NbQ11g09350.1 NbQ11g13660.1 NbQ11g13930.1 NbQ11g16260.1 NbQ11g17360.1 NbQ11g18430.1 NbQ11g21080.1 NbQ11g23280.1 NbQ11g23990.1 NbQ11g25050.1 NbQ12g03770.1 NbQ12g04850.1 NbQ12g07340.1 NbQ12g09080.1 NbQ12g10820.1 NbQ12g12070.1 NbQ12g14750.1 NbQ12g15000.1 NbQ12g15230.1 NbQ12g20380.1 NbQ12g21080.1 NbQ12g21830.1 NbQ12g23960.1 NbQ13g01300.1 NbQ13g02350.1 NbQ13g03860.1 NbQ13g04410.1 NbQ13g08800.1 NbQ13g09850.1 NbQ13g10370.1 NbQ13g11700.1 NbQ13g12420.1 NbQ13g15780.1 NbQ13g16040.1 NbQ13g23160.1 NbQ13g24120.1 NbQ13g24540.1 NbQ13g25080.1 NbQ13g25490.1 NbQ13g28240.1 NbQ13g29770.1 NbQ14g01070.1 NbQ14g03950.1 NbQ14g05360.1 NbQ14g05410.1 NbQ14g06880.1 NbQ14g07270.1 NbQ14g07500.1 NbQ14g10290.1 NbQ14g10770.1 NbQ14g14320.1 NbQ14g17890.1 NbQ14g18710.1 NbQ14g20960.1 NbQ14g22890.1 NbQ15g00150.1 NbQ15g02300.1 NbQ15g02330.1 NbQ15g02350.1 NbQ15g03230.1 NbQ15g06190.1 NbQ15g07120.1 NbQ15g07750.1 NbQ15g09000.1 NbQ15g09050.1 NbQ15g11920.1 NbQ15g12650.1 NbQ15g12840.1 NbQ15g15670.1 NbQ15g15930.1 NbQ15g18670.1 NbQ15g19070.1 NbQ15g20620.1 NbQ15g22880.1 NbQ15g23000.1 NbQ15g26060.1 NbQ16g00880.1 NbQ16g04360.1 NbQ16g06490.1 NbQ16g09100.1 NbQ16g11020.1 NbQ16g11560.1 NbQ16g13810.1 NbQ16g13820.1 NbQ16g17040.1 NbQ16g17130.1 NbQ16g17340.1 NbQ16g18390.1 NbQ16g18430.1 NbQ16g23100.1 NbQ16g23570.1 NbQ16g24270.1 NbQ16g25200.1 NbQ16g25830.1 NbQ16g25880.1 NbQ16g25990.1 NbQ16g26610.1 NbQ16g26660.1 NbQ16g28010.1 NbQ16g28180.1 NbQ17g01150.1 NbQ17g01180.1 NbQ17g01570.1 NbQ17g01950.1 NbQ17g05460.1 NbQ17g05540.1 NbQ17g05980.1 NbQ17g07990.1 NbQ17g08300.1 NbQ17g09330.1 NbQ17g09400.1 NbQ17g10090.1 NbQ17g11220.1 NbQ17g13030.1 NbQ17g15460.1 NbQ17g16690.1 NbQ17g20980.1 NbQ17g22370.1 NbQ17g25040.1 NbQ17g28730.1 NbQ18g02140.1 NbQ18g02740.1 NbQ18g05440.1 NbQ18g06120.1 NbQ18g07470.1 NbQ18g12320.1 NbQ18g12530.1 NbQ18g12850.1 NbQ18g13840.1 NbQ18g14420.1 NbQ18g14930.1 NbQ18g15730.1 NbQ18g17750.1 NbQ18g17850.1 NbQ18g21060.1 NbQ19g01040.1 NbQ19g05480.1 NbQ19g06450.1 NbQ19g06510.1 NbQ19g08330.1 NbQ19g11840.1 NbQ19g11880.1 NbQ19g13750.1 NbQ19g14190.1 NbQ19g14210.1 NbQ19g14920.1 NbQ19g18540.1 NbQ19g19870.1 NbQ19g21020.1 NbQ19g21220.1 NbQ19g22080.1 NbQ19g22800.1 NbQ19g24690.1 NbQ19g24730.1 rna19561
OG0000001 Capann_59V1aChr01g048170.1 NbL01g00020.1 NbL01g00940.1 NbL01g02330.1 NbL01g03550.1 NbL01g03650.1 NbL01g04410.1 NbL01g04920.1 NbL01g06850.1 NbL01g16120.1 NbL01g19150.1 NbL01g20140.1 NbL01g20930.1 NbL01g22230.1 NbL01g24190.1 NbL01g24280.1 NbL01g24300.1 NbL02g00570.1 NbL02g00900.1 NbL02g01270.1 NbL02g02110.1 NbL02g02210.1 NbL02g02470.1 NbL02g03180.1 NbL02g04740.1 NbL02g04750.1 NbL02g06120.1 NbL02g06860.1 NbL02g07280.1 NbL02g07680.1 NbL02g07740.1 NbL02g09780.1 NbL02g11320.1 NbL02g12670.1 NbL02g13080.1 NbL02g14050.1 NbL02g14190.1 NbL02g15010.1 NbL02g15890.1 NbL02g16190.1 NbL02g16730.1 NbL02g17070.1 NbL02g17360.1 NbL02g18820.1 NbL02g19340.1 NbL02g20100.1 NbL02g23950.1 NbL02g24800.1 NbL03g01610.1 NbL03g01680.1 NbL03g01890.1 NbL03g02230.1 NbL03g02600.1 NbL03g03410.1 NbL03g04990.1 NbL03g05400.1 NbL03g08030.1 NbL03g08250.1 NbL03g08690.1 NbL03g10230.1 NbL03g11060.1 NbL03g13030.1 NbL03g14960.1 NbL03g15110.1 NbL03g16690.1 NbL03g16900.1 NbL03g18260.1 NbL03g18950.1 NbL03g21180.1 NbL03g21210.1 NbL03g21530.1 NbL03g22960.1 NbL03g24430.1 NbL04g01140.1 NbL04g01490.1 NbL04g02030.1 NbL04g02560.1 NbL04g03700.1 NbL04g04160.1 NbL04g05240.1 NbL04g05420.1 NbL04g05850.1 NbL04g12420.1 NbL04g12640.1 NbL04g13650.1 NbL04g13780.1 NbL04g14310.1 NbL04g16260.1 NbL04g17750.1 NbL04g18380.1 NbL04g18870.1 NbL04g19030.1 NbL04g19630.1 NbL05g00320.1 NbL05g03060.1 NbL05g03300.1 NbL05g04060.1 NbL05g07620.1 NbL05g08630.1 NbL05g09580.1 NbL05g10060.1 NbL05g11400.1 NbL05g12280.1 NbL05g13170.1 NbL05g16020.1 NbL05g17530.1 NbL05g17730.1 NbL05g18340.1 NbL05g18590.1 NbL05g18600.1 NbL05g19640.1 NbL05g20640.1 NbL05g21000.1 NbL05g21640.1 NbL05g22610.1 NbL06g00640.1 NbL06g00660.1 NbL06g02210.1 NbL06g03150.1 NbL06g03680.1 NbL06g04910.1 NbL06g07950.1 NbL06g09970.1 NbL06g11480.1 NbL06g12220.1 NbL06g12400.1 NbL06g12460.1 NbL06g12850.1 NbL06g13120.1 NbL06g14450.1 NbL06g14780.1 NbL06g16990.1 NbL06g17200.1 NbL06g17760.1 NbL06g20380.1 NbL07g02100.1 NbL07g02540.1 NbL07g02970.1 NbL07g03110.1 NbL07g04840.1 NbL07g05350.1 NbL07g06580.1 NbL07g07530.1 NbL07g08450.1 NbL07g09380.1 NbL07g09870.1 NbL07g10730.1 NbL07g10850.1 NbL07g11080.1 NbL07g12450.1 NbL07g12710.1 NbL07g13110.1 NbL07g13920.1 NbL07g14240.1 NbL07g15520.1 NbL07g16220.1 NbL07g17480.1 NbL08g01820.1 NbL08g02750.1 NbL08g02930.1 NbL08g03510.1 NbL08g03620.1 NbL08g03850.1 NbL08g03970.1 NbL08g04040.1 NbL08g06150.1 NbL08g06410.1 NbL08g06680.1 NbL08g06730.1 NbL08g07620.1 NbL08g08450.1 NbL08g08640.1 NbL08g09910.1 NbL08g10160.1 NbL08g11760.1 NbL08g12570.1 NbL08g13630.1 NbL08g13890.1 NbL08g15050.1 NbL08g15340.1 NbL08g18010.1 NbL08g18420.1 NbL08g19080.1 NbL08g19190.1 NbL09g00900.1 NbL09g02160.1 NbL09g02330.1 NbL09g02470.1 NbL09g04150.1 NbL09g05210.1 NbL09g07010.1 NbL09g09070.1 NbL09g10290.1 NbL09g10500.1 NbL09g11220.1 NbL09g13490.1 NbL09g15290.1 NbL09g15830.1 NbL09g17240.1 NbL09g19250.1 NbL09g19460.1 NbL09g20190.1 NbL09g21040.1 NbL09g21520.1 NbL09g23460.1 NbL10g00080.1 NbL10g03710.1 NbL10g04330.1 NbL10g04560.1 NbL10g05200.1 NbL10g06320.1 NbL10g07510.1 NbL10g07960.1 NbL10g08670.1 NbL10g08970.1 NbL10g11120.1 NbL10g11340.1 NbL10g11820.1 NbL10g13720.1 NbL10g14560.1 NbL10g14770.1 NbL10g16430.1 NbL10g18140.1 NbL10g18380.1 NbL10g19280.1 NbL10g19690.1 NbL10g21210.1 NbL10g22680.1 NbL10g23160.1 NbL10g23560.1 NbL10g24210.1 NbL11g00680.1 NbL11g00970.1 NbL11g01230.1 NbL11g01270.1 NbL11g01520.1 NbL11g01530.1 NbL11g02920.1 NbL11g03540.1 NbL11g03990.1 NbL11g05630.1 NbL11g08950.1 NbL11g08980.1 NbL11g09510.1 NbL11g10840.1 NbL11g11030.1 NbL11g11230.1 NbL11g12430.1 NbL11g13300.1 NbL11g15430.1 NbL11g16390.1 NbL11g16410.1 NbL11g17320.1 NbL11g18090.1 NbL11g21310.1 NbL11g21470.1 NbL11g21780.1 NbL11g21820.1 NbL11g22270.1 NbL11g22310.1 NbL11g23180.1 NbL11g24100.1 NbL12g00130.1 NbL12g01810.1 NbL12g02230.1 NbL12g02720.1 NbL12g02760.1 NbL12g04120.1 NbL12g04550.1 NbL12g06630.1 NbL12g07830.1 NbL12g09170.1 NbL12g10580.1 NbL12g12090.1 NbL12g12490.1 NbL12g12630.1 NbL12g12800.1 NbL12g13320.1 NbL12g13460.1 NbL12g14430.1 NbL12g14970.1 NbL12g15490.1 NbL12g17460.1 NbL12g18190.1 NbL12g18590.1 NbL12g19900.1 NbL12g20690.1 NbL12g22040.1 NbL12g22560.1 NbL13g00350.1 NbL13g01440.1 NbL13g02400.1 NbL13g03210.1 NbL13g03360.1 NbL13g04070.1 NbL13g05250.1 NbL13g08460.1 NbL13g09010.1 NbL13g09140.1 NbL13g10290.1 NbL13g11570.1 NbL13g13370.1 NbL13g14910.1 NbL13g18680.1 NbL13g19510.1 NbL13g23520.1 NbL13g24010.1 NbL13g24190.1 NbL13g24460.1 NbL13g26310.1 NbL13g26640.1 NbL13g26860.1 NbL13g27260.1 NbL13g27960.1 NbL14g02460.1 NbL14g02750.1 NbL14g08750.1 NbL14g08910.1 NbL14g09120.1 NbL14g09540.1 NbL14g09920.1 NbL14g11070.1 NbL14g11150.1 NbL14g12570.1 NbL14g14530.1 NbL14g14860.1 NbL14g15240.1 NbL14g15460.1 NbL14g16620.1 NbL14g16910.1 NbL14g17940.1 NbL14g21150.1 NbL14g21750.1 NbL14g21910.1 NbL15g00790.1 NbL15g01170.1 NbL15g02310.1 NbL15g04220.1 NbL15g05970.1 NbL15g06340.1 NbL15g06440.1 NbL15g07020.1 NbL15g07370.1 NbL15g07470.1 NbL15g09010.1 NbL15g13210.1 NbL15g14550.1 NbL15g14600.1 NbL15g17290.1 NbL15g18170.1 NbL15g19710.1 NbL15g21840.1 NbL15g21930.1 NbL15g23410.1 NbL15g23420.1 NbL15g23430.1 NbL15g25130.1 NbL15g25200.1 NbL16g01760.1 NbL16g02140.1 NbL16g04460.1 NbL16g05010.1 NbL16g05020.1 NbL16g06780.1 NbL16g07540.1 NbL16g07980.1 NbL16g09760.1 NbL16g10610.1 NbL16g12320.1 NbL16g13510.1 NbL16g14420.1 NbL16g15690.1 NbL16g17420.1 NbL16g17790.1 NbL16g17880.1 NbL16g18730.1 NbL16g18940.1 NbL16g19440.1 NbL16g20980.1 NbL16g23180.1 NbL16g23610.1 NbL16g23660.1 NbL16g23910.1 NbL16g24550.1 NbL16g24640.1 NbL16g25300.1 NbL16g25630.1 NbL16g26710.1 NbL17g01590.1 NbL17g02070.1 NbL17g02120.1 NbL17g02920.1 NbL17g03040.1 NbL17g03540.1 NbL17g03700.1 NbL17g03800.1 NbL17g05400.1 NbL17g07510.1 NbL17g08450.1 NbL17g08930.1 NbL17g10090.1 NbL17g14370.1 NbL17g14600.1 NbL17g15390.1 NbL17g15900.1 NbL17g16000.1 NbL17g16910.1 NbL17g17480.1 NbL17g18240.1 NbL17g20020.1 NbL17g20830.1 NbL17g21220.1 NbL17g21690.1 NbL17g25960.1 NbL18g00030.1 NbL18g00150.1 NbL18g00310.1 NbL18g00670.1 NbL18g00700.1 NbL18g01630.1 NbL18g02650.1 NbL18g04460.1 NbL18g05210.1 NbL18g05690.1 NbL18g07270.1 NbL18g07440.1 NbL18g07500.1 NbL18g09090.1 NbL18g09810.1 NbL18g09880.1 NbL18g10500.1 NbL18g10990.1 NbL18g12070.1 NbL18g13060.1 NbL18g17480.1 NbL19g00230.1 NbL19g04470.1 NbL19g04700.1 NbL19g07770.1 NbL19g07870.1 NbL19g09260.1 NbL19g09870.1 NbL19g13480.1 NbL19g14360.1 NbL19g14400.1 NbL19g14720.1 NbL19g17700.1 NbL19g18860.1 NbL19g22750.1 NbQ13g00370.1 NbQ14g10310.1 NbQ17g06050.1

This script appears to have some problems with the file:
import pandas as pd

orthofinder_output = "../OrthoFinder-res/Results_Jan23/Orthogroups/Orthogroups-fixed.txt"
orthogroups = pd.read_csv(orthofinder_output, sep=' ', header=None)

# Extract gene family information
expansions = {}
contractions = {}
for i, row in orthogroups.iterrows():
    if len(row) > 2:
        expansions[row[0]] = row[2:]
    else:
        contractions[row[0]] = row[1]

# Print expansions and contractions
print("Expansions:", expansions)
print("Contractions:", contractions)

I got the following error:
geneFamilyExpansionsContractions.py:
 20: DtypeWarning:
   Columns (2,3,4,5,6,...,964,965,966,967,968) have mixed types.
   Specify dtype option on import or set low_memory=False.
  orthogroups = pd.read_csv(orthofinder_output, sep=' ', header=None)
...
Expansions:
...
Name: 63241, Length: 967, dtype: object, 'OG0063242':
2      NaN
3      NaN
4      NaN
5      NaN
6      NaN
      ...
964    NaN
965    NaN
966    NaN
967    NaN
968    NaN
Name: 63242, Length: 967, dtype: object}
Contractions: {}

How can it be fixed?

The answer is actually right there in the error message (warning message, actually):

DtypeWarning: Columns (2,... 968) have mixed types. Specify dtype option on import or set low_memory=False.

Setting low_memory-False makes the error go away:
orthogroups = pd.read_csv(orthofinder_output, sep=' ', header=None, low_memory=False)

According to this random blog post I found when I searched for "pandas DtypeWarning: have mixed types", this is because:

Pandas is really nice because instead of stopping altogether, it guesses which dtype a column has. However, by default, Pandas has the low_memory=True argument. This means that the CSV file gets split up in multiple chunks and the guess is done for every chunk, resulting in a column with multiple dtypes. Pandas is so kind to let you know it was confused and something might have happened.

The same blog also suggests that this isn't really the right approach, and that makes sense, since this is forcing Pandas to read the entire file and that might be a problem for huge files or low memory machines. So, reading on, the blog suggests the other option provided by the error, setting the dtype manually.
Now, I am afraid I have never used Pandas and don't really know much Python either, but a quick glance at the Pandas documentation suggests that the object data type is the most versatile and should work for everything, including strings. Since your file is nothing but a collection of strings, as far as I can tell, this seems like the best approach.
So I changed your script to this:
#!/usr/bin/env python3

import pandas as pd
import sys

orthofinder_output = sys.argv[1]
orthogroups = pd.read_csv(orthofinder_output, sep=' ', header=None, low_memory=False)
# Extract gene family information
expansions = {}
contractions = {}
for i, row in orthogroups.iterrows():
    if len(row) > 2:
        expansions[row[0]] = row[2:]
    else:
        contractions[row[0]] = row[1]

# Print expansions and contractions
print("Expansions:", expansions)
print("Contractions:", contractions)

And it seems to work just fine with no warnings shown.

The error right at the end is

Specify dtype option on import or set low_memory=False.

Edit Oh wait @terdon has already spotted this! :-D ... great minds and all that :)
In addition, the pandas print out it states ...
dtype object

This isn't good. Name: 63242, is a bit unusual BTW worth checking.
Looking at the values ...
NaN

... that is really not good.
What I think is happening is it's been assigned to a numerical datatype, but it's categorical data. Hence NaN because like the alphabet can be none-base 10 numbers, in data compression used a lot, but it's got to be consistent pattern (I genuinely don't know if pandas can spot/work in none-base 10 number systems). pandas is really C (probably C++) - its a static type backend - so those datatypes are the key to this.
It makes sense because your data is categorical data and NaN means its managed to confuse a large amount of it for e.g. int64.
Looking the documentation the alternative is
orthogroups = pd.read_csv(orthofinder_output, sep=' ', header=None, dtype='category')

The documentation is

low_memorybool, default True
Internally process the file in chunks, resulting in lower memory use while parsing, but possibly mixed type inference. To ensure no mixed types either set False, or specify the type with the dtype parameter. Note that the entire file is read into a single DataFrame regardless, use the chunksize or iterator parameter to return the data in chunks. (Only valid with C parser).

What you want is to issue a df.dtypes command once it is loaded and the output will define the datatype per column,
'A'    int64
'B'    int32
'C'   category

If you have mixed dtypes it may (I stress may) be able to correct it via,
df.astype('int64') # this is not your data
df.astype('category') # this is your data

HOWEVER, I don't think this will work, because NaN indicates its trashed the data. Anyway the data is definitely categorical BTW.
Anyway, good question and interesting stuff.

As a bit of gossip, you can also switch to categorical ordered data I think this will work for dataframes, this example is Series. There's other ways to order the data BTW.
from pandas.api.types import CategoricalDtype
my_dtype = CategoricalDtype(\
    categories=[2, 1], ordered=True)
df.dtype(my_dtype)

The output is 1 and 2 now switch places because 2 < 1. Categorial datatypes look to be the right datatype for your inputs.


text:
Remove low quality reads
I want to remove reads from FASTQ file that contain homopolymers > 10bp and remove reads with <35 average quality score across the entire read.
To remove homopolymers > 10bp, I tried this on a Linux machine, but it only removes the sequence line:
zcat file.fastq.gz |
  awk '!/A{10,}/&&!/C{10,}/&&!/G{10,}/&&!/T{10,}/ {print}' cleaned_file.fastq

Please [edit] your question and explain what you need in more detail. What do you mean by "less than the 35 quality score"? Do you want to remove all reads that have even a single base with <35 score? Do you want to calculate the average quality across the entire read and remove reads with <35 average quality score? Something else? Why do you want to remove these from the fastq file instead of handling it in the aligned bam or as part of your variant calling process? How are you defining homopolymers here? Do you just want to remove reads if they have >10 consecutive identical bases?

Use an off the shelf tool for read preprocessing.
Here is one:
./fastq_qual_trimmer -i test.fq -m 35 -H 10

That one does exactly what you want but seems a little old/unmaintained, so here is another that does more or less the same thing:
fastq-mcf --qual-mean 35 --homopolymer-pct {X} adapters.fa reads.fq

where {X} is 10 / read length, adapters.fa is an adapter file (which I believe can be empty or filled with dummy sequences).
You could also use a library like biopython or dnaio to write a quick script to do this, but it hardly seems worth it.

Thanks for your help. Sorry I didn't write in detail but i need to calculate the average quality over the whole read and remove the reads with average quality score <35. As I understand it, these tools can't handle it.

@n_n sorry, used wrong options. Both of these tools can handle it across the read. Updated answer.

Thank you so much for taking the time to help me. Your assistance is greatly appreciated.

@user_01 if this answer has solved your problem could you very kindly please mark this as accepted? It would be a small way of showing gratitude and helps the site in several ways.

Just one more question I tried fastq-mcf as you suggest for my fastq files but i got this Error: number of input files must match number of '-o' output files. May you help me to understand whats wrong with code.
   read_length=$(head -n 2 ${file} | tail -n 1 | wc -m)
   homopolymer_pct=$((8/${read_length}))
   fastq-mcf ${file} adapters.fa --qual-mean 35 --homopolymer-pct  ${homopolymer_pct} -o filtered/${file}
done```

@user_01 I don't have the time right now to extensively troubleshoot but i'd guess that something is wrong with the way that you're using {$file}, or with the files you are passing, or with the existence of directory filtered/. If you print out {$file} or some of the other values you might see something.


text:
Does the contrast matrix I have made address the question I am trying to answer?
I have the following design matrix:
mm_noreps.interactions <- model.matrix(~condition*TRAPed)

Both variables are factors condition has 4 levels and TRAPed has 2 levels.
pheatmap(mm_noreps.interactions)

I want to use the following contrasts:

TRAP positive (just consider this a drug treatment) vs TRAP negative in all conditions,
c(-1,0,0,0,0,1,0,0)
FearRecallpositive vs all other condition's positive cells c(0,0,-1/3,0,-1/3,-1/3,0,1), and
FearRecallpostive vs all other conditions c(-1/7,-1/7,-1/7,-1/7,-1/7,-1/7,-1/7,1).

The main question I have is are there differentially expressed genes specific to the FearRecall:TRAPedpostive interaction.  Does contrast 3 address this? The other two are more sanity checks than actual contrasts of interest.