text:
Error while running computeMatrix command in Deeptools
I am trying to get computeMatrix for bigwig file in specific genomic region using deeptools.
Below is the code I am using
computeMatrix reference-point --referencePoint TSS \
-b 1000 -a 1000 \
-R ~/Desktop/ATAC/ATAC/Inducible_elements_Greenberg.bed \
-S ~/Desktop/ATAC/Control.mRp.clN.bigWig \
--skipZeros \
-o ~/Desktop/ATAC/matrix_controlmRPATACBasalcomputematrix.gz \

But I get the following error
File "/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/numpy/__init__.py", line 305, in __getattr__
    raise AttributeError(__former_attrs__[attr])
AttributeError: module 'numpy' has no attribute 'float'.
`np.float` was a deprecated alias for the builtin `float`. To avoid this error in existing code, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.
The aliases was originally deprecated in NumPy 1.20; for more details and guidance see the original release note at:
    https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations. Did you mean: 'cfloat'?

These are the packages in my Conda environment
# packages in environment at /Users/abhinavsoni/anaconda3/envs/ATAC_Seq:
#
# Name                    Version                   Build  Channel
bzip2                     1.0.8                h3422bc3_4    conda-forge
ca-certificates           2023.5.7             hf0a4a13_0    conda-forge
libblas                   3.9.0           16_osxarm64_openblas    conda-forge
libcblas                  3.9.0           16_osxarm64_openblas    conda-forge
libcxx                    16.0.3               h4653b0c_0    conda-forge
libffi                    3.4.2                h3422bc3_5    conda-forge
libgfortran               5.0.0           12_2_0_hd922786_31    conda-forge
libgfortran5              12.2.0              h0eea778_31    conda-forge
liblapack                 3.9.0           16_osxarm64_openblas    conda-forge
libopenblas               0.3.21          openmp_hc731615_3    conda-forge
libsqlite                 3.41.2               hb31c410_1    conda-forge
libzlib                   1.2.13               h03a7124_4    conda-forge
llvm-openmp               16.0.3               h1c12783_0    conda-forge
ncurses                   6.3                  h07bb92c_1    conda-forge
numpy                     1.23.0          py310h0a343b5_0    conda-forge
openssl                   3.1.0                h53f4e23_3    conda-forge
pip                       23.1.2             pyhd8ed1ab_0    conda-forge
python                    3.10.11         h3ba56d0_0_cpython    conda-forge
python_abi                3.10                    3_cp310    conda-forge
readline                  8.2                  h92ec313_1    conda-forge
setuptools                67.7.2             pyhd8ed1ab_0    conda-forge
tk                        8.6.12               he1e0b03_0    conda-forge
tzdata                    2023c                h71feb2d_0    conda-forge
wheel                     0.40.0             pyhd8ed1ab_0    conda-forge
xz                        5.2.6                h57fd34a_0    conda-forge

What could be its possible work around?

It worked, I installed a the newer version of deeptools, from GitHub:
https://github.com/deeptools/deepTools


text:
CNV Kit ` from . import commands ImportError: cannot import name 'commands' from '__main__'`
I am trying to run some code for my colleague in bash
/path/to/cnvkit.py batch /path/to/my/folder/with/bams/*.bam \
  --normal --targets ${bed_file} \
  --fasta path/to/my/resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta \
  --output-reference /path/to/my/CD_BATCH1_reference.cnn \
  --output-dir /path/to/my/Group_1 --scatter

however, I keep getting this truly peculiar error
Traceback (most recent call last):
  File "/path/to/cnvkit.py", line 4, in <module>
    from . import commands
ImportError: cannot import name 'commands' from '__main__' (/path/to/cnvkit.py)

Attempt 1
I checked if the path was correct and multiple times I found that it was correct. It did not change the error.
Attempt 2
Perhaps I need to chmod +x * everything in the CNV kit directory and that will help. It did not change the error.
Why is this error occurring? How do I solve this issue?

Could you download it from conda and try again? Its here, https://docs.conda.io/en/main/miniconda.html
Then,
conda create --name myenv
conda activate myenv
conda install -c bioconda cnvkit
cnvkit batch /path/to/my/folder/with/bams/*.bam --normal --targets ${bed_file} --fasta path/to/my/resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta --output-reference /path/to/my/CD_BATCH1_reference.cnn --output-dir /path/to/my/Group_1 --scatter

The path to the file is not used, the command is simply cnvkit
This avoids all the privilege stuff and ensures you're running a virtual environment to protect your system Python.

To address the comments. "conda" is slow and doesn't ever load.
There is an alternative to conda called mamba, this is fast and designed to address this exact issue. I have never needed to use it so the following is a rough guide only.
Download and install via
conda create --name myenv
conda activate myenv
conda install -c conda-forge mamba

The documentation is here, https://mamba.readthedocs.io/en/latest/user_guide/mamba.html
I've never used it but the idea is ...
mamba create -n mambaenv <list of packages>
mamba activate mambaenv

The list of packages could be left blank for now, then,
mamba install "Python>=3.10.0" numpy pandas pathlib
mamba install cnvkit # SEE NOTE

Can also do
mamba install -n mambaenv cnvkit

NOTE cnvkit exists strictly in -c bioconda. You do not seem to search and install via channels in mamba
You search for a package and its dependencies via
mamba repoquery search cnvkit
mamba repoquery depends -a cnvkit

It is important to place a virtual environment between the system Python and package Python (don't want a system corruption). More to the question, this would solve all the permission problems without changing your system in anyway (kinda important).

I did try this before and I tried it again now but for whatever reason I get Solving environment: failed with initial frozen solve. Retrying with flexible solve. over and over again. It never works though for the command conda install -c bioconda cnvkit.

Okay @Indira thanks. Addressed above

"I have never needed to use it" you lucky so-and-so! We regular joes need to stick to mamba if we need conda environments solved in one human lifetime.

@RamRS this is funny. I suggest the reason is that "Regular Joses" (pronounced Hosays) where I am, don't use conda. "Regular Joes" where you are, all use conda. Thus your server always times out. :-)

source: Stack Exchange

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

sort: votes


text:
NCBI Blast databases - cannot find some databases
I'm working with a system that does local BLAST searches using some preconfigured BLAST databases.  From what I see, these were updated within the last year, but I cannot find any current databases available that have those same names.  One specific set is a non-redundant set of human protein or nucleotide sequences.  The databases are named human_nr and human_nt (with various extensions as per the BLAST db format).  There are others names up_sp_human (Uniprot-Swissprot human?) and up_tbl_human (Uniprot-Tremble human?), but a cannot find any similar databases on the NCBI site or other sites.
Where else should I look for such preconfigured databases?

I don't think there is a human_nr or human_nt database for NCBI.
An equivalent database is accessible via the standard nr/nt database and simply restrict the organism parameter to human or to taxa id 9606 . This essentially will do the same thing.
Or maybe this sort of thing?,
blastdbcmd -db nr -entry all -outfmt "%g %T" | awk ' { if ($2 == 9606) { print $1 } } ' > human_flat-list.db

This is super helpful!  I'm guessing these were custom prepared similar to what you describe.  Thank you!

Thanks for your feedback @KirkD-CO

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


text:
Find Lipidic Mediator receptors
I want to study the ligand-receptor interaction and I am searching for lipidic ligands( mediators), I have found 16 so far. Does anyone know any database where I could find more lipidic ligands. I was searching on LIPID MAPS website but I couldn't know how to filter search for the lipidic ligands. Any recommendation is appreciated.


text:
Assign score to PyRanges from another PyRanges object
I have two PyRanges objects. PR1 is a set of regions of interest across a chromosome, PR2 tiles the genome and has a score associated with each range.
I need a way to get for each range in PR1 its score based on its overlap with PR2. If it is contained within just one range in PR2, it inherits the score of that range. If it overlaps multiple ranges on PR2, the score should be a weighted average of those ranges based on the size of the overlap.

Applying this to the image above, it would be:
score(A) = score(X)
score(B) = (OverlapBX * score(X) + OverlapBY * score(Y)) / sizeB

I've been browsing the PyRanges docs and I see a lot of functions that do almost that, but most of them don't take into account a score.
I had done this before in R using GRanges, where I could iterate over the overlap object to piece together scores, but I can't find a direct port of R's findOverlaps() in PyRanges.

The inputs I'm using for this are a list of mutation positions and a map of recombination rates over the chromosome they're on. So if two mutations are in a section with a constant recombination rate, the section between them inherits that recombination rate, but if the map says there's different rates between them, I need to take that into account.
Here's a pastebin with the data I'm working with: https://pastebin.com/sBYsjRQB

The Unfun Cat's first suggestion:
dist_to_pos = distances.join(positions, report_overlap=True).drop(like="_b")
df = dist_to_pos.df
g = df.groupby("Id2")
avg_score = g.cMperMb.sum() / g.Overlap.first()
avg_score.name = "AvgScore"
j_with_avg_score = pr.PyRanges(df.set_index("Id2").join(avg_score))

The results this give are... curious. Here's what it gives when applied to my data:
j_with_avg_score

Start   cMperMb     cMatstart   End     Chromosome  Id1     Overlap     AvgScore
0   2916415     1.223681    3.272236    2917484     chr1    0   1   12.236813
1   2916415     1.223681    3.272236    2917484     chr1    1   24  12.236813
2   2916415     1.223681    3.272236    2917484     chr1    2   5   12.236813
3   2916415     1.223681    3.272236    2917484     chr1    3   74  12.236813
4   2916415     1.223681    3.272236    2917484     chr1    4   110     12.236813
...     ...     ...     ...     ...     ...     ...     ...     ...
947     2960324     0.269208    3.291457    2960462     chr1    913     138     0.001951
948     2960462     0.352644    3.291494    2960790     chr1    913     77  0.018319
949     2960462     0.352644    3.291494    2960790     chr1    914     6   0.018319
950     2960462     0.352644    3.291494    2960790     chr1    915     15  0.018319
951     2960462     0.352644    3.291494    2960790     chr1    916     71  0.018319

952 rows × 8 columns
## input data contains 917 sections --> mismatch!

res_dic = {}
for x,y in zip(j_with_avg_score.Start, j_with_avg_score.AvgScore):
    res_dic[x] = y
res_dic

{2916415: 12.236813262999997,
 2917484: 0.02265139573888889,
 2917538: 0.06125764894,
 2917608: 0.10651243958260867,
 2918281: 0.6792060563708334,
 2918836: 1.1527867390640016,
 2921187: 4.176119918333334,
 2921274: 1.2543028034,
 2921300: 0.9058800047555555,
 2921976: 3.26162101122,
 2922306: 0.4522015679581395,
 2923593: 0.7409518310923077,
 2924172: 8.386433179800006,
 2925614: 0.31748233525000014,
 2926340: 0.17712928483599996,
 2926675: 1.4567026153333333,
 2926730: 5.063891720272731,
 2927172: 2.42991683355,
 2927410: 0.9370814839,
 2927411: 0.1303501307571428,
 2928125: 0.20534225654799998,
 2928557: 0.19228154158888885,
 2928849: 1.5796778519999994,
 2929269: 1.9946231915199983,
 2931597: 0.15826380113170732,
 2933043: 0.010025524,
 2934387: 0.05559167711694925,
 2938888: 0.0023231045055555556,
 2940194: 0.01882597932,
 2941104: 0.1670264778,
 2941694: 0.0197562495886076,
 2942700: 0.3931709166833331,
 2946526: 1.6788888500999994,
 2947460: 0.19597836537931007,
 2951262: 5.532616244500001,
 2956899: 0.44889279189999975,
 2960324: 0.0019507828376811595,
 2960462: 0.018319161953246753}

Something's getting chewed up here, due to the number of unique sections in the output, the duplicates, and the start and end positions not matching the input data.
However, the first value calculated is exactly the same as what I had gotten in my R script (with a few magnitude adjustments): 2916415: 12.236813262999997 vs 2917178: 0.00012236813263
So the calculations happening are most likely correct, I think it's mainly the start and end positions that are getting messed up.

I think join will get you halfway there? It is meant to be the replacement for findOverlaps

@TheUnfunCat I had started to fiddle with joins but then tried to use cluster and merge afterwards, but I got some numpy errors from those. I'll log them on the github.

@TheUnfunCat tried to replicate the error but it works now. I can't remember the details but it was something about numpy not recognising the variable type 'long'

Have you tried updating PyRanges?
Anyways, you can do gr.df to get the dataframe afterward your join and then do your calculations in pandas.

I have some inspiration below. It computes a weighted score based on overlap.
I need to know what exactly you want the overlap BX and BY to be. Could you post some example data and expected output?
# initialize test data
import numpy as np
np.random.seed(0)

import pyranges as pr
a1 = pr.data.aorta().drop()
a1.Id1 = list(range(len(a1)))
a2 = pr.data.aorta2().drop()
a2.Id2 = list(range(len(a2)))
a1.Score1 = np.random.random(len(a1)) * 100
a2.Score2 = np.random.random(len(a2)) * 100

# compute weighted average based on overlaps
j = a1.join(a2, report_overlap=True).drop(like="_b")

df = j.df

g = df.groupby("Id2")
avg_score = g.Score2.sum() / g.Overlap.first()
avg_score.name = "AvgScore"
j_with_avg_score = pr.PyRanges(df.set_index("Id2").join(avg_score))

Ideally I would like to combine BX and BY to B(XY), and get the average score during that process. I'll update the question with the files I'm using.

I tried the method, but it didn't work great. I'll update the post with the info from that


text:
Refactoring pandas using an iterator via chunksize
This question was also asked on Stack Overflow
Bioinformatics rationale eggNOG files can be very big and sump all available RAM for regular to medium sized desktops.

I am looking for advice on using a pandas iterator.
I performed a parsing operation using Python pandas, the size of the input files is  eggNOG a program used in metagenomic gene function identification but has broader application. The files are large and parsing is resulting in 'RAM bottleneck' phenomenon and will not processing/parse the file.
The obvious solution is to shift to an iterator, which for pandas is the  chunksize option
import pandas as pd
import numpy as np

df = pd.read_csv('myinfile.csv', sep="\t", chunksize=100)

Whats changed with the original code is the chunksize=100 bit, forcing an iterator.
The next step is just to perform a simple operation, dropping a few columns and moving all '-' characters to np.nan then writing the whole file to disk.
df.drop(['score', 'evalue', 'Description', 'EC', 'PFAMs'],axis=1).replace('-', np.nan)
df.to_csv('my.csv',sep='\t',index=False)

How is the step performed under an iterator?

The solution posted by @el_oso and @Corralien here was as follows:

if you're already having memory issues reading the file, the below will actually read and write to the new file in chunks.

cols_to_keep = {'col1', 'col2', 'col3'}
df = pd.read_csv('myinfile.csv', chunksize=100, usecols=cols_to_keep)
with open('my.csv', 'w') as f:
    f.write('\t'.join(cols_to_keep) + '\n')
    for chunk in df:
        #do your processing here, appending that chunk to the file
        chunk.to_csv(f, header=None, sep='\t')

f is _io.TextIOWrapper i.e. a stream you can consistently write to in small memory chunks

This is very elegant because there are two components which are both cool and highlight why pandas is so widely used.

Instead of loading a large file and then dumping stuff thats not wanted: don't load junk at point source.
Simply open the write file outside the chunking loop will then stream the chunks to straight to disk and once the final chunking loop completes the file snaps shut.

The process is

Load a chunk of data, but don't load any 'junk'.

Process the chunk of data and write it to disk.

Load another chunk of data repeating the above two steps.

Continue until the entire file is processed.

The "do processing" is an OO module (there's alot going on) and the chunk is actually self.chunk. None of the original code needs to change its simply restructured inside the open loop to stream chunks to disk.
The beauty of this solution is chunksize is a variable, here 100 is just the example. Python's iterator called yield is one line at a time iteration, but this the automates a defined 'chunk' at a time. Thus the runtime speed can be tuned against the "RAM bottleneck".
Testing The chunked outfile was then reduced via groupby to remove the "duplicates" and frequencies summed (code not shown). The chunking OO module output file is identical to "non-chunked" OO module and written back into the module.
Its worth upvoting the SO answers.


text:
change the order in ggplot
I'm trying to make a plot that looks like this:

This is a rudimentary plot I did in excel. It's a plot of reads by length. Very simple nothing crazy to it. The important thing to note is that it's plotting the data in increasing order of length.
When I try to use ggplot to make this same plot in R it does something like this:

Where it orients the read names in alphabetical order not in ascending order of the lengths. The dataframe/table I have read into R is already listed in ascending length order. For example:
1   c6dd66b4-dab5-484b-92b6-c4ecdcc155b1   1570
2   ecf15d4f-15b1-4636-a97d-e0928281e311   1739
3   461f05c7-0fc6-4564-95df-6188966fe616   1790
4   e989c32c-491e-4af1-a37d-4d4f5c87a3d9   1797
5   ce94c2b3-442b-4b8c-bd59-6d7c23d5fe3d   1824
6   ed15ea16-a28d-42ca-9a07-2818d90ddd4b   1832
7   481550cd-6f87-44c7-915f-f0f3d638ffe4   1853

As you can see the lengths increase as you go down and the read names are in no particular order. Can someone help me read in the correct order to make this plot actually work in ggplot?

the technical term you want is called "ranking". The question (and the answer) will be much clearer once that is used. rank is one of the possible R methods that will do this.

You need to enforce the order you need on the read-name column. Let's say df is your data frame with read.name and length as the two columns.
This is a fairly easy operation.
df <- df[order(df$length),]
df$read.name <- factor(df$read.name, levels = df$read.name, ordered = TRUE)
ggplot(df, ...) #Your current ggplot code goes here

Hmm this gave me an error factor level [99] is duplicated but I guess I should just google that too, right? No offense but there's no need to be rude with your response. I did Google this question extensively before making this post and the order and factor methods did pop up and I did try them but didn't get it to work which is why I posted this....

My apologies @rimo I'm a bit slow on this. The solution to your data is that it needs ranking before its plotted, for example df$rank <-  rank(df$length) then if you place the df$rank into ggplot2 that should work. RamRS is using order instead.  If it was Python I'd code it for you (coding is easier than writing). It is fairly easy with R however.

Yes, you should google that too. That's how we learn. Google error messages and try stuff until you either understand the whole picture or run into something a lot beyond your understanding. The error says your read names are duplicated. In that case, I'd either make them unique by adding a new component, use rank as M__ suggests or plot length vs 1:nrow(df) in the x axis so sorting by x-axis doesn't do a thing, as your plot doesn't really need read names on an axis.


text:
h5ad file format filter
I'm trying to do simple filter for single data that is stored in h5ad file format using this
import scanpy as sc

sc.pp.filter_cells(adata, min_genes=500)
sc.pp.filter_genes(adata, min_cells=10)
print(adata.n_obs, adata.n_vars)

when I saved the filtered data and see the size of my original data compare it, the size of the filtered data is more, but the dimension etc are reduced due to filtering.
So is it something wrong am I doing?
UPDATE
adata = sc.read_h5ad('/Single_cell/cerebellar_development/508cf892-174c-45ab-b2dc-05f54f1ee7ed/aldinger20.processed.h5ad')
adata

AnnData object with n_obs × n_vars = 69174 × 49166
    obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'sample_id', 'percent.mito', 'S.Score', 'G2M.Score', 'Phase', 'CC.Difference', 'nCount_SCT', 'nFeature_SCT', 'age', 'figure_clusters', 'sex', 'type', 'experiment', 'fig_cell_type'
    var: 'name'
    obsm: 'X_pca', 'X_umap', 'X_tsne'

sc.pp.filter_cells(adata, min_genes=500)
sc.pp.filter_genes(adata, min_cells=10)
print(adata.n_obs, adata.n_vars)

results_file = 'adata.h5ad'
adata.write(results_file)

FULL CODE
updated version: h5ad filter, added layers to make the filtered file usable for seurat to read and make a correct seurat object
#https://github.com/satijalab/seurat/issues/2485

import scanpy as sc
import numpy as np
import pandas as pd
import time
from ipywidgets import interact, interactive, widgets
from IPython.display import display, HTML

def filter_and_save_h5ad(h5ad_file, min_genes, min_cells):
    # h5ad file read
    adata = sc.read_h5ad(h5ad_file)

    # cells and genes filter parameter from user
    sc.pp.filter_cells(adata, min_genes=min_genes)
    sc.pp.filter_genes(adata, min_cells=min_cells)

    # Add 'counts' layer
    adata.layers['counts'] = adata.X.copy()

    # add 'winsorization' layer
    adata.layers['winsorized'] = adata.X.copy()

    # filtered dataset save
    filtered_file = h5ad_file.replace('.h5ad', '_filtered.h5ad')
    adata.write(filtered_file)

    print(f"Filtered dataset saved to {filtered_file}")

# input widgets
h5ad_file_widget = widgets.Text(description='File_add:')
min_genes_widget = widgets.IntText(description='Min_Genes:')
min_cells_widget = widgets.IntText(description='Min_Cells:')
button = widgets.Button(description='Filter_Save')

# timer widget
time_widget = widgets.HTML(value="
Time passed away: 0 seconds
")

def on_button_clicked(click_me):
    h5ad_file = h5ad_file_widget.value.strip()
    min_genes = min_genes_widget.value
    min_cells = min_cells_widget.value

    # timer
    start_time = time.time()

    filter_and_save_h5ad(h5ad_file, min_genes, min_cells)

    # Calculate and update the elapsed time
    elapsed_time = time.time() - start_time
    time_widget.value = f"
Time passed away: {elapsed_time} seconds
"

button.on_click(on_button_clicked)

# Display the widgets and button
display(h5ad_file_widget, min_genes_widget, min_cells_widget, button, time_widget)

You didn't show us how the original .h5ad object was generated, how you read it in, how you wrote it, and if you did anything else in-between besides filtering. AnnData.write() takes a compression arg that might explain the difference. As also mentioned in that link, ensuring you are writing a sparse matrix can make a substantial difference to file size (the matrix in adata.X is generally by far the largest part of the object)

I was taking this dataset https://data.humancellatlas.org/explore/projects/85a9263b-0887-48ed-ab1a-ddfa773727b6/get-curl-command

I have updated the question

I have updated my code please have a look and for this "You didn't show us how the original .h5ad object was generated" please have a look https://data.humancellatlas.org/explore/projects/85a9263b-0887-48ed-ab1a-ddfa773727b6

Scanpy's h5ad is in fact hdf5 format. This format is used for rapid read/write (i/o) large amounts of data and thats its purpose. Its fairly easy to use from scratch.
Its purpose is not to compress data. If you are concerned with the contents of h5ad why not just read them directly via the hdf5 viewer. You can download it here. https://www.hdfgroup.org/downloads/hdfview/
The file will then open up as a "Explorer" (Windows) or "Finder" (OSX) GUI file system and you can open the various files and see whats in them.
My suspicion is that the original data is been written into the file. 'cause it can i/o rapidly and access precisely what it needs then file size could be a secondary consideration. Anyway you can peer in and find out by looking at the original file and the hdf5 save you've done.
An example from the website is here

Cool right?
Personally, I think HDF5 is overrated and feather outperforms it for speed (i/o) - but thats just for dataframes. Its not much faster than pickling in case you are interested and pickling in addition does compression. parquet is the format that does rapid i/o but specialises in compression, again mostly for dataframes.
HDF5 however does permit very pretty GUI interfaces.

Also if you want to load up and look at various files, that just uses the Python hdf5 module.


text:
Plotting infercnv results
I'm working with matched single cell data, where we have treated and untreated samples for the same patient. I ran CNV analysis using the infercnv package.
I've followed the tutorial:
# data matrix
counts_matrix <- scData@assays$RNA@counts
meta = data.frame(labels = Idents(scData), row.names = names(Idents(scData)))
unique(meta$labels) # check the cell labels

# Create the infercnv object
infercnv_obj <- CreateInfercnvObject(raw_counts_matrix = counts_matrix,
                                     annotations_file = meta,
                                     delim = "\t",
                                     gene_order_file = "hg38_gencode_v27.txt",
                                     ref_group_names = NULL)

# perform infercnv operations to reveal cnv signal
infercnv_obj = infercnv::run(infercnv_obj,
                             cutoff=1,  # use 1 for smart-seq, 0.1 for 10x-genomics
                             out_dir="output_dir",  # dir is auto-created for storing outputs
                             cluster_by_groups=T,   # cluster
                             denoise=T,
                             HMM=T
                             )

Now, I'm not sure how to plot the results. Any help would be appreciated.
When running the below recommendation I get this:
INFO [2023-05-19 12:49:22] ::plot_cnv:Start
INFO [2023-05-19 12:49:22] ::plot_cnv:Current data dimensions (r,c)=34108,17624 Total=82463956 Min=0 Max=55099.
INFO [2023-05-19 12:49:22] ::plot_cnv:Depending on the size of the matrix this may take a moment.
<sparse>[ <logic> ]: .M.sub.i.logical() maybe inefficient
Warning: sparse->dense coercion: allocating vector of size 4.5 GiBINFO [2023-05-19 12:50:12] plot_cnv(): auto thresholding at: (-2.725632 , 3.000000)
INFO [2023-05-19 12:50:35] plot_cnv_observation:Start
INFO [2023-05-19 12:50:35] Observation data size: Cells= 17624 Genes= 34108
INFO [2023-05-19 12:50:35] clustering observations via method: ward.D
INFO [2023-05-19 12:50:35] Number of cells in group(1) is 331
INFO [2023-05-19 12:50:36] group size being clustered:  331,34108
Error in parallelDist(data_to_cluster, threads = infercnv.env$GLOBAL_NUM_THREADS) :
  x must be a matrix or a list of matrices.

Personally I would not use intercnv_obj throughout, because is there was a bug in line 2, i.e. using = instead of <- then the created intercnv_obj might be trying to be passed into plot_cnv.
Could you try
data(infercnv_obj)

plot_cnv(
  infercnv_obj,
  out_dir = ".",
  title = "inferCNV",
  obs_title = "Observations (Cells)",
  ref_title = "References (Cells)",
  cluster_by_groups = TRUE,
  cluster_references = TRUE,
  plot_chr_scale = FALSE,
  chr_lengths = NULL,
  k_obs_groups = 1,
  contig_cex = 1,
  x.center = mean(infercnv_obj@expr.data),
  x.range = "auto",
  hclust_method = "ward.D",
  custom_color_pal = NULL,
  color_safe_pal = FALSE,
  output_filename = "infercnv",
  output_format = "png",
  png_res = 300,
  dynamic_resize = 0,
  ref_contig = NULL,
  write_expr_matrix = FALSE,
  write_phylo = FALSE,
  useRaster = TRUE
)

If the above doesn't work ...
What I think is happening is the input data is the problem. This precise error is already known and was reported on Broad's GitHub here:
https://github.com/broadinstitute/infercnv/issues/415
The responder (Broad Institute) believed it was the version problem, i.e. the version was outdated. However a separate respondee believed it was the input data format and suggested shifting away from a sparse matrix. That would make sense in context to the error being thrown.

Hi , I get an error when running that. I updated my question with the message.

@mmpp I have modified the code and I've included an error log on GitHub. Apologies for the delay. If problem continued I'd recommend referring to this GitHub error log ... https://github.com/broadinstitute/infercnv/issues/415 or else filing a separate error and referring to this error log. It is not clear whether Broad has fully resolved this - but they already know about it. Good luck :-)

What did you change?

data(infercnv_obj) prior plot_cnv

It’s absolutely fine to use = as an assignment operator in R, it’s largely a myth that <- is ‘better’ in some way. It’s certainly not a ‘bug’ by any means.

If you haven't already, I recommend checking out the infercnv YouTube tutorial (https://github.com/broadinstitute/inferCNV/wiki/Running-InferCNV) for detailed instructions on running InferCNV. Additionally, you can explore Uphyloplot2 (https://github.com/harbourlab/uphyloplot2) for building a tree from the output generated by InferCNV.
I noticed that you didn't use the untreated sample as the normal reference for your CNV analysis. I'm curious to understand your reasoning behind this. In my case, I also have an untreated sample that isn't diploid, along with two treated samples. Initially, I planned to use the untreated sample as the "normal" reference. However, considering that InferCNV assumes a normal diploid reference, I'm now uncertain if using the untreated control would be appropriate.
Please let me know your thoughts, any help would be much appreciated!


text:
How to associate point mutations to quantitative protein levels
My experiment is based on a quantitative continuous variable levels (B12 vitamin) in a sample of patients retrospectively selected based on values that can take (under our criteria):

normal
low levels: under the physiological or normal reference values but not enough to be considered deficiency
deficiency (which is associated with poor nutrition, gastric disease, etc.). These data are excluded.

We are interested in low levels. Our hypothesis is to prove if mutations in 10 genes can be associated with these physiological low levels. These physiological low levels are quantified under the format of maximum and minimum value of B12. The genes can take 3 possible categories:

wt: absence of mutation
het: half of mutation
hom: full mutation

Hypothesis:

I would like, if possible, to try something like a genetic risk score, if multiple het or hom states in multiple genes can be additive. I don't know how to approach this, because there can be multiple permutations, and maybe permutation 1 has a stronger association with low levels than permutation 2.
Any ideas on how to carry this out? I don't really know if mutations are additive or opposite.
Not sure if some clustering approach can be helpful here, like PCA or euclidean distance, linear discriminant analysis... Any advice is welcome.

My idea is to execute a summary of levels of B12 for each gene and see in a boxplot possible differences, perhaps look for some association (or correlation), but I doubt it.
Any idea or previous input or experiment like this would be tremendously welcome.

The solution is supervised machine learning (ML) and is fairly standard analytics for precisely this problem.  PCA or euclidean distance, will not prove this especially distance methods. Discriminant analysis is okay. What you need are complex regression analytics and that broadly includes ML. So IMO  the project needs an analyst. I hope thats helpful and again its just my opinion

Could you recommend me something to begin?

I personally would recommend collaboration, unless someone already published an algorithm specific for purpose. Its a bit like a specialist/expert in military logistics wanting to perform a "commando assault", i.e. it would be really hard, particularly if the key source of help was online.

Yes, you can run an additive model (i.e., B*AA = 2*B*Aa, aa = B*AA - 4*Aa) or a dominance model (B*AA = B*Aa). You can test both of these and see which is a better fit to your data.
However, it is practically known already that the model will be at least predominantly additive.


text:
How to normalise indel of different transcript direction?
I am trying to reproduce an analysis result, and there is a requirement regarding the format of the result: the '3' rule ': all mutations should be expressed in a position close to the end of the gene transcription direction (3' end).
The result I obtained is(annotate from VEP)
Chromosome  Start_Position  End_Position    Strand  Variant_Classification  Variant_Type    Reference_Allele    Tumor_Seq_Allele1   Tumor_Seq_Allele2   HGVSc   HGVSp   HGVSp_Short Transcript_ID
10  89692835    89692840    +   In_Frame_Del    DEL GATCTT  GATCTT  -   c.321_326del    p.Leu108_Asp109del  p.L108_D109del  NM_000314.8

the standard answer
chr10   89692837    89692842    TCTTGA      PTEN    NM_00314    c.321_326del    p.Leu108_Asp109del

I speculate that ：
the HGVS standard followed by VEP includes the 3'rule . So ultimately, there is the same p.Leu108_ Asp109del。 But the output of vep did not adjust Start_ Position and End_ Position？（Perhaps it thinks that after obtaining the correct amino acid mutation , it doesn't matter anymore）
But for me, I hope the ultimate Start_ Position and End_ The Position also complies with the regulations, which seems a bit troublesome because the different  transcripts may have different directions and always rely on the IGV to correct it , which is also tiring.
Thanks for any suggestions

I have added an IGV image to better illustrate the problem. Reference genome regions 89692830 to 89692845 for the example in the figure.
bwa(or other software，like bowtie2) believes that the position of del should be
ref:GTGAAGATCTTGACCA
alt:GTGAA------GACCA

But 3'rule believes that the position of del should be
ref:GTGAAGATCTTGACCA
alt:GTGAAGA------CCA

To verify whether the results match 3'rule, it is necessary to know the direction of the transcripts, and each gene always has multiple transcripts, which is where I am confused

Thanks for the update @tanqiang

Can you clarify which representation you need? Is it the genomic level, the coding level or the protein level? Also, what is the input here? Do you start with variants or with VEP output? How many do you have? Do you only want to get one result, or does this need to be automated? Why are other transcripts relevant if you already have NM_000314.8?


text:
How to get a GISAID account? I registered months ago, still no reply!
Inspired by amateur variant hunters, I would like to join the Pango lineage proposal community and help contribute to variant surveillance.
However, I cannot seem to get access to GISAID, the platform that exclusively hosts a large share of SARS-CoV-2 sequences from around the world.
I applied/registered for an account back in summer 2022 (almost a year ago) using this form: https://gisaid.org/register/
But since then, nothing has happened. I have send an email and filled out the contact form on the website but never heard back.
I am very confused how GISAID can say they are a "public" database, while at the same time not giving accounts to people like me, not even responding to emails.
How does everyone else get access? Is there a secret email you need to contact? Or do you need to pay? I've read here which describes the founder of GISAID in detail.
Any help would be greatly appreciated!

I've got a GISAID account and it's subject to an approval process. It was pretty simple and is awarded quickly. There are definitely no fees needed. I personally didn't see the advantages over NCBI's really very impressive SARS-CoV-2 database system and thats what I'd recommend.
Their resource is here:
https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/virus?SeqType_s=Nucleotide&VirusLineage_ss=SARS-CoV-2,%20taxid:2697049&ids=MZ039713%20MZ039682%20OV522182%20
Its perfect for interrogation, it has the full spectrum of epidemiological metadata, and further details are in this post here:
How to programatically download SARS-CoV-2 fasta from NCBI/Genbank via API?
The Science journal page seemed a very good write up of GISAID BTW.

If you are wanting GISAID for networking, thats a bit tricky because NCBI is not going to do that.
You can reapply for GISAID access, but at the time I was much more impressed by NCBI. The speed at which NCBI produced an awesome database was truly jaw-dropping and they'd implemented full pangolineages classification at an early stage. Thus I just use NCBI.

Thanks for the edit @gringer, much appreciated. ".. it has all the data" is good interpretation of what I wrote, but my writing wasn't clear. Specifically, "it has all the epidemiological metadata" is what I intended to say.


text:
Define and calculate pep_sc and pep_sc_noref in pyrosetta
I can calculate total rosetta energy score of a peptide-protein complex. However, when using flexpepdock module from pyrosetta I didn't find anything direct parameter like pep_sc or pep_sc_noref in the pyrosetta module. I try to write a workaround for pep_sc.
from pyrosetta import rosetta

# Initialize PyRosetta
pyrosetta.init()

# Load the protein-peptide complex from the PDB file
pose = rosetta.core.import_pose.pose_from_file("./all_pdb/clean_pdb/1d4t.clean.pdb")

# Create a score function
scorefxn = rosetta.core.scoring.ScoreFunction()

# Set up the FlexPepDock protocol
flexpep_dock = rosetta.protocols.flexpep_docking.FlexPepDockingProtocol()

# Apply the FlexPepDock protocol
flexpep_dock.apply(pose)

# Get the energy terms of interest
score = pose.energies().total_energy()
hbond_sc = pose.energies().residue_total_energies(1)[rosetta.core.scoring.hbond_sc]
rama = pose.energies().residue_total_energies(1)[rosetta.core.scoring.rama]
pep_sc = pose.energies().residue_total_energies(2)[rosetta.core.scoring.fa_atr] + \
         pose.energies().residue_total_energies(2)[rosetta.core.scoring.fa_rep] + \
         pose.energies().residue_total_energies(2)[rosetta.core.scoring.fa_sol] + \
         pose.energies().residue_total_energies(2)[rosetta.core.scoring.fa_elec]

fa_atr = pose.energies().residue_total_energies(1)[rosetta.core.scoring.fa_atr]
fa_pair = pose.energies().residue_total_energies(1)[rosetta.core.scoring.fa_pair]
#pep_sc_noref = pose.energies().residue_total_energies(1)[rosetta.core.scoring.pep_sc_noref]
fa_rep = pose.energies().residue_total_energies(1)[rosetta.core.scoring.fa_rep]
fa_dun = pose.energies().residue_total_energies(1)[rosetta.core.scoring.fa_dun]
fa_sol = pose.energies().residue_total_energies(1)[rosetta.core.scoring.fa_sol]
#startRMSbb = pose.energies().residue_total_energies(1)[rosetta.core.scoring.startRMSbb]
# Calculate I_sc (Interface energy)
complex_energy = pose.energies().total_energies()[rosetta.core.scoring.total_score]
protein_energy = pose.energies().residue_total_energies(1)[rosetta.core.scoring.total_score]
peptide_energy = pose.energies().residue_total_energies(2)[rosetta.core.scoring.total_score]
I_sc = complex_energy - (protein_energy + peptide_energy)

# Print the energy terms
print("Score:", score)
print("hbond_sc:", hbond_sc)
print("rama:", rama)
print("I_sc:", I_sc)
print("pep_sc:", pep_sc)
print("fa_atr:", fa_atr)
print("fa_pair:", fa_pair)
print("fa_rep:", fa_rep)
print("fa_dun:", fa_dun)
print("fa_sol:", fa_sol)

Can anyone tell me how to calculate pep_sc and pep_sc_noref for this code? Also, can anyone tell me what is the formula for pep_sc and pep_sc_noref if there is no direct reference .


text:
How to use Gblocks for trimming single-copy gene sequences?
I have run OrthoFinder on a set of 11 genomes and got the results in a folder. From the output folder, I saved all the single-copy gene sequences in a single FASTA file and now wish to remove all gaps using Gblocks before using them for constructing a phylogenomic tree. But I don't know the command line used for running Gblocks. Any help will be appreciated!

I assume you've aligned these sequences? I don't know GBlocks but I know trimal and will do the same thing (or do it better). Could you describe what you mean by 'removing gaps'? 1. Are you looking to trim the 5' and 3' ends? 2. Wanting to remove 'gap [indel] only columns'? 3. May be wanting to remove the any position which contains an indel, i.e. the entire site is removed? The last option is the best approach.

Did you look for the gblocks docs? You know it's >20 years old that tool!

Hey @M__ ! I am still learning the procedures for doing phylogenetic analysis but I think I need to remove ambiguous sequences that hamper the accuracy of a phylogeny. You can check this paper which has cited Gblocks.                                        [link] (https://doi.org/10.1016/j.syapm.2020.126166)          If you feel trimal is better could you please explain how to use it? Thanks!

@Chris_Rands, I did come across the Gblocks documentation but it did not work for me. I think it has not been updated in a while. Since my supervisor shared the paper which used Gblocks, I wanted to learn it. I agree that the latest tool would be better.

Basically this is 16S data. If its proteins there's a different strategy - even if the data set needed is nucleotides (separate question). I would simply use trimal which is available on conda
trimal -in <inputfile> -out <outputfile> -nogaps

The input file is aligned fasta here and will output as aligned fasta.
Any site (alignment position) with an indel - will be removed. This is definitely preferred if performing an iqtree analysis, because of the way it handles indels. This will also clean up poorly aligned, alignments.
e.g.
---ATG
---A-G
--TATG

This will result in...
AG
AG
AG

Ultimately for alignments - there's no magic bullet - and manual editing and inspection is expected e.g. using seqotron. There are very good reasons for this concerning circular mathematics. However, modern aligners such as muscle5 are extremely good.
Trimal is available via conda
conda install -c bioconda trimal

In terms of using GBlock I'd need to know whether it was simply deleting indel sites, deleting indel sites at >50% frequency or whether it was performing a genetic diversity calculation per site. I don't think it is doing the latter. It is better to perform diversity saturation calculation explicitly rather than using an unknown automate because thats kinda important. If its deleting sites with excessive >50% indels - honestly its better to do a blockwise deletion rather than do something half-half. Within SSU there is a reason might be being done because there can be excessive indels at deep taxonomic levels with this data type. Again using iqtree just do it 'blockwise' which is the trimal output.
Note I've not read the original GBlock paper, but if you summarise its application I could consider it against trimal. However, I definitely prefer trimal considering the downstream application.

Thank you again, @M__ . Will definitely go ahead with what you have suggested. Although, I must admit, being an absolute amateur in bioinformatics I don't follow some of the things you have so sincerely explained. Forgive me if I am being imprudent, do you have any way that I could learn all this under your guidance? I wish to be thorough.. thanks again!

If you need the theory you can always ask a separate question. The theory is very detailed, why certain things and done and otherwise not.


text:
shiny app for single cell cluster export issue
I'm trying to create a shiny app for single cell data visualization and then label cells with genes till this step it works, the final step when I try to export cluster it is not exporting the desired cluster or i might be doing something wrong which Im not sure
This the code
library(Seurat)
library(BiocParallel)
library(reticulate)
library(Matrix)
library(shiny)
library(dplyr)

# Set the maximum file upload size to 3 GB
options(shiny.maxRequestSize = 3e+9)

# Define the UI
ui <- fluidPage(
  titlePanel("Single-cell RNA-seq Analysis"),
  sidebarLayout(
    sidebarPanel(
      fileInput("h5ad_file", "Upload H5AD File", accept = ".h5ad"),
      actionButton("run_button", "Run Analysis"),
      selectInput("gene_input", "Select Gene for Feature Plot:", choices = NULL),
      downloadButton("export_button", "Export Cluster")
    ),
    mainPanel(
      plotOutput("feature_plot")
    )
  )
)

# Define the server
server <- function(input, output, session) {

  # Function to perform the analysis
  performAnalysis <- function(h5ad_file) {
    sc <- import("scanpy")
    adata <- sc$read_h5ad(h5ad_file)

    counts <- t(adata$layers["counts"])
    colnames(counts) <- adata$obs_names$to_list()
    rownames(counts) <- adata$var_names$to_list()
    counts <- Matrix::Matrix(as.matrix(counts), sparse = TRUE)

    data <- t(adata$layers["winsorized"])
    colnames(data) <- adata$obs_names$to_list()
    rownames(data) <- adata$var_names$to_list()
    data <- Matrix::Matrix(as.matrix(data), sparse = TRUE)

    seurat <- CreateSeuratObject(counts)
    seurat <- SetAssayData(seurat, "data", data)
    seurat <- AddMetaData(seurat, adata$obs)

    seurat[["percent.mt"]] <- PercentageFeatureSet(seurat, pattern = "^MT[-\\.]")
    seurat <- NormalizeData(seurat)
    seurat <- FindVariableFeatures(seurat, nfeatures = 500)
    seurat <- ScaleData(seurat)
    seurat <- ScaleData(seurat, vars.to.regress = c("nFeature_RNA", "percent.mt"))
    seurat <- RunPCA(seurat, npcs = 50)
    seurat <- RunUMAP(seurat, dims = 1:20)
    seurat <- FindNeighbors(seurat, dims = 1:20)
    seurat <- FindClusters(seurat, resolution = 1)
    cl_markers <- FindAllMarkers(seurat, only.pos = TRUE, min.pct = 0.25, logfc.threshold = log(1.2))
    top_gene <- cl_markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)

    return(seurat)
  }

  # Reactive value to store the Seurat object
  seurat_object <- reactiveVal(NULL)

  # Perform analysis when the Run Analysis button is clicked
  observeEvent(input$run_button, {
    req(input$h5ad_file)
    h5ad_file <- input$h5ad_file$datapath
    seurat <- performAnalysis(h5ad_file)
    seurat_object(seurat)

    # Update gene choices for feature plot
    gene_choices <- rownames(seurat@assays$RNA@counts)
    updateSelectInput(session, "gene_input", choices = gene_choices)
  })

  # Render the feature plot
  output$feature_plot <- renderPlot({
    seurat <- seurat_object()
    if (!is.null(seurat)) {
      gene <- input$gene_input
      FeaturePlot(seurat, features = gene, label = TRUE)
    }
  })

  # Export cluster when the Export Cluster button is clicked
  observeEvent(input$export_button, {
    seurat <- seurat_object()
    if (!is.null(seurat)) {
      cluster_data <- cbind(Cluster = seurat@meta.data$cluster, seurat@assays$RNA@counts)
      write.csv(cluster_data, "cluster.csv", row.names = TRUE)
      updateActionButton(session, "export_button", label = "Download Cluster", href = "cluster.csv")
    }
  })

  # Function to end the session
  endSession <- function() {
    stopApp()  # This will end the Shiny session and shut down the server
  }

  # Call the endSession function when the session ends (browser window closed)
  session$onSessionEnded(endSession)
}

# Run the Shiny app
shinyApp(ui = ui, server = server)

The seurat object which is generated its structure is like this
str(seurat)
Formal class 'Seurat' [package "SeuratObject"] with 13 slots
  ..@ assays      :List of 1
  .. ..$ RNA:Formal class 'Assay' [package "SeuratObject"] with 8 slots
  .. .. .. ..@ counts       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. .. .. ..@ i       : int [1:2238732] 29 73 80 148 163 184 186 227 229 230 ...
  .. .. .. .. .. ..@ p       : int [1:2639] 0 779 2131 3260 4220 4741 5522 6304 7094 7626 ...
  .. .. .. .. .. ..@ Dim     : int [1:2] 13714 2638
  .. .. .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ...
  .. .. .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. .. .. ..@ x       : num [1:2238732] 1 1 2 1 1 1 1 41 1 1 ...
  .. .. .. .. .. ..@ factors : list()
  .. .. .. ..@ data         :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. .. .. ..@ i       : int [1:2238732] 29 73 80 148 163 184 186 227 229 230 ...
  .. .. .. .. .. ..@ p       : int [1:2639] 0 779 2131 3260 4220 4741 5522 6304 7094 7626 ...
  .. .. .. .. .. ..@ Dim     : int [1:2] 13714 2638
  .. .. .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ...
  .. .. .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. .. .. ..@ x       : num [1:2238732] 1.64 1.64 2.23 1.64 1.64 ...
  .. .. .. .. .. ..@ factors : list()
  .. .. .. ..@ scale.data   : num [1:500, 1:2638] 1.5265 -0.4982 -0.0324 -0.3129 -0.4631 ...
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : chr [1:500] "MRPL20" "RER1" "TNFRSF9" "SRM" ...
  .. .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. ..@ key          : chr "rna_"
  .. .. .. ..@ assay.orig   : NULL
  .. .. .. ..@ var.features : chr [1:500] "PPBP" "LYZ" "S100A9" "IGLL5" ...
  .. .. .. ..@ meta.features:'data.frame':  13714 obs. of  5 variables:
  .. .. .. .. ..$ vst.mean                 : num [1:13714] 0.00341 0.00114 0.0019 0.00114 0.00682 ...
  .. .. .. .. ..$ vst.variance             : num [1:13714] 0.0034 0.00114 0.00189 0.00114 0.00678 ...
  .. .. .. .. ..$ vst.variance.expected    : num [1:13714] 0.00365 0.00114 0.00197 0.00114 0.00748 ...
  .. .. .. .. ..$ vst.variance.standardized: num [1:13714] 0.933 0.992 0.963 0.992 0.906 ...
  .. .. .. .. ..$ vst.variable             : logi [1:13714] FALSE FALSE FALSE FALSE FALSE FALSE ...
  .. .. .. ..@ misc         : list()
  ..@ meta.data   :'data.frame':    2638 obs. of  10 variables:
  .. ..$ orig.ident     : Factor w/ 1 level "SeuratProject": 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ nCount_RNA     : num [1:2638] 2419 4903 3147 2639 980 ...
  .. ..$ nFeature_RNA   : int [1:2638] 779 1352 1129 960 521 781 782 790 532 550 ...
  .. ..$ n_genes        : num [1:2638] 781 1352 1131 960 522 ...
  .. ..$ percent_mito   : num [1:2638] 0.0302 0.0379 0.0089 0.0174 0.0122 ...
  .. ..$ n_counts       : num [1:2638] 2419 4903 3147 2639 980 ...
  .. ..$ louvain        : Factor w/ 8 levels "CD4 T cells",..: 1 3 1 2 5 4 4 4 1 6 ...
  .. ..$ percent.mt     : num [1:2638] 3.02 3.79 0.89 1.74 1.22 ...
  .. ..$ RNA_snn_res.1  : Factor w/ 8 levels "0","1","2","3",..: 1 3 1 6 4 1 5 5 5 6 ...
  .. ..$ seurat_clusters: Factor w/ 8 levels "0","1","2","3",..: 1 3 1 6 4 1 5 5 5 6 ...
  ..@ active.assay: chr "RNA"
  ..@ active.ident: Factor w/ 8 levels "0","1","2","3",..: 1 3 1 6 4 1 5 5 5 6 ...
  .. ..- attr(*, "names")= chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  ..@ graphs      :List of 2
  .. ..$ RNA_nn :Formal class 'Graph' [package "SeuratObject"] with 7 slots
  .. .. .. ..@ assay.used: chr "RNA"
  .. .. .. ..@ i         : int [1:52760] 0 102 187 274 457 833 1147 1503 1506 1550 ...
  .. .. .. ..@ p         : int [1:2639] 0 17 24 30 33 40 61 83 85 96 ...
  .. .. .. ..@ Dim       : int [1:2] 2638 2638
  .. .. .. ..@ Dimnames  :List of 2
  .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. ..@ x         : num [1:52760] 1 1 1 1 1 1 1 1 1 1 ...
  .. .. .. ..@ factors   : list()
  .. ..$ RNA_snn:Formal class 'Graph' [package "SeuratObject"] with 7 slots
  .. .. .. ..@ assay.used: chr "RNA"
  .. .. .. ..@ i         : int [1:245992] 0 69 102 172 187 189 248 274 276 315 ...
  .. .. .. ..@ p         : int [1:2639] 0 79 118 179 253 323 424 541 574 650 ...
  .. .. .. ..@ Dim       : int [1:2] 2638 2638
  .. .. .. ..@ Dimnames  :List of 2
  .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. ..@ x         : num [1:245992] 1 0.0811 0.2903 0.1765 0.1765 ...
  .. .. .. ..@ factors   : list()
  ..@ neighbors   : list()
  ..@ reductions  :List of 2
  .. ..$ pca :Formal class 'DimReduc' [package "SeuratObject"] with 9 slots
  .. .. .. ..@ cell.embeddings           : num [1:2638, 1:50] -3.544 -5.892 -3.544 10.139 0.932 ...
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. .. .. ..$ : chr [1:50] "PC_1" "PC_2" "PC_3" "PC_4" ...
  .. .. .. ..@ feature.loadings          : num [1:500, 1:50] 0.02701 0.16566 0.17077 -0.00786 -0.03674 ...
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : chr [1:500] "PPBP" "LYZ" "S100A9" "IGLL5" ...
  .. .. .. .. .. ..$ : chr [1:50] "PC_1" "PC_2" "PC_3" "PC_4" ...
  .. .. .. ..@ feature.loadings.projected: num[0 , 0 ]
  .. .. .. ..@ assay.used                : chr "RNA"
  .. .. .. ..@ global                    : logi FALSE
  .. .. .. ..@ stdev                     : num [1:50] 5 3.39 3.2 2.78 1.94 ...
  .. .. .. ..@ key                       : chr "PC_"
  .. .. .. ..@ jackstraw                 :Formal class 'JackStrawData' [package "SeuratObject"] with 4 slots
  .. .. .. .. .. ..@ empirical.p.values     : num[0 , 0 ]
  .. .. .. .. .. ..@ fake.reduction.scores  : num[0 , 0 ]
  .. .. .. .. .. ..@ empirical.p.values.full: num[0 , 0 ]
  .. .. .. .. .. ..@ overall.p.values       : num[0 , 0 ]
  .. .. .. ..@ misc                      :List of 1
  .. .. .. .. ..$ total.variance: num 457
  .. ..$ umap:Formal class 'DimReduc' [package "SeuratObject"] with 9 slots
  .. .. .. ..@ cell.embeddings           : num [1:2638, 1:2] -3.44 -7.94 -1.77 8.79 -1.51 ...
  .. .. .. .. ..- attr(*, "scaled:center")= num [1:2] 1.64 -1.3
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. .. .. ..$ : chr [1:2] "UMAP_1" "UMAP_2"
  .. .. .. ..@ feature.loadings          : num[0 , 0 ]
  .. .. .. ..@ feature.loadings.projected: num[0 , 0 ]
  .. .. .. ..@ assay.used                : chr "RNA"
  .. .. .. ..@ global                    : logi TRUE
  .. .. .. ..@ stdev                     : num(0)
  .. .. .. ..@ key                       : chr "UMAP_"
  .. .. .. ..@ jackstraw                 :Formal class 'JackStrawData' [package "SeuratObject"] with 4 slots
  .. .. .. .. .. ..@ empirical.p.values     : num[0 , 0 ]
  .. .. .. .. .. ..@ fake.reduction.scores  : num[0 , 0 ]
  .. .. .. .. .. ..@ empirical.p.values.full: num[0 , 0 ]
  .. .. .. .. .. ..@ overall.p.values       : num[0 , 0 ]
  .. .. .. ..@ misc                      : list()
  ..@ images      : list()
  ..@ project.name: chr "SeuratProject"
  ..@ misc        : list()
  ..@ version     :Classes 'package_version', 'numeric_version'  hidden list of 1
  .. ..$ : int [1:3] 4 1 3
  ..@ commands    :List of 7
  .. ..$ NormalizeData.RNA       :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
  .. .. .. ..@ name       : chr "NormalizeData.RNA"
  .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2023-05-21 03:58:21"
  .. .. .. ..@ assay.used : chr "RNA"
  .. .. .. ..@ call.string: chr "NormalizeData(seurat)"
  .. .. .. ..@ params     :List of 5
  .. .. .. .. ..$ assay               : chr "RNA"
  .. .. .. .. ..$ normalization.method: chr "LogNormalize"
  .. .. .. .. ..$ scale.factor        : num 10000
  .. .. .. .. ..$ margin              : num 1
  .. .. .. .. ..$ verbose             : logi TRUE
  .. ..$ FindVariableFeatures.RNA:Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
  .. .. .. ..@ name       : chr "FindVariableFeatures.RNA"
  .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2023-05-21 03:58:22"
  .. .. .. ..@ assay.used : chr "RNA"
  .. .. .. ..@ call.string: chr "FindVariableFeatures(seurat, nfeatures = 500)"
  .. .. .. ..@ params     :List of 12
  .. .. .. .. ..$ assay              : chr "RNA"
  .. .. .. .. ..$ selection.method   : chr "vst"
  .. .. .. .. ..$ loess.span         : num 0.3
  .. .. .. .. ..$ clip.max           : chr "auto"
  .. .. .. .. ..$ mean.function      :function (mat, display_progress)
  .. .. .. .. ..$ dispersion.function:function (mat, display_progress)
  .. .. .. .. ..$ num.bin            : num 20
  .. .. .. .. ..$ binning.method     : chr "equal_width"
  .. .. .. .. ..$ nfeatures          : num 500
  .. .. .. .. ..$ mean.cutoff        : num [1:2] 0.1 8
  .. .. .. .. ..$ dispersion.cutoff  : num [1:2] 1 Inf
  .. .. .. .. ..$ verbose            : logi TRUE
  .. ..$ ScaleData.RNA           :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
  .. .. .. ..@ name       : chr "ScaleData.RNA"
  .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2023-05-21 03:58:23"
  .. .. .. ..@ assay.used : chr "RNA"
  .. .. .. ..@ call.string: chr [1:2] "ScaleData(seurat, vars.to.regress = c(\"nFeature_RNA\", " "    \"percent.mt\"))"
  .. .. .. ..@ params     :List of 11
  .. .. .. .. ..$ features          : chr [1:500] "PPBP" "LYZ" "S100A9" "IGLL5" ...
  .. .. .. .. ..$ assay             : chr "RNA"
  .. .. .. .. ..$ vars.to.regress   : chr [1:2] "nFeature_RNA" "percent.mt"
  .. .. .. .. ..$ model.use         : chr "linear"
  .. .. .. .. ..$ use.umi           : logi FALSE
  .. .. .. .. ..$ do.scale          : logi TRUE
  .. .. .. .. ..$ do.center         : logi TRUE
  .. .. .. .. ..$ scale.max         : num 10
  .. .. .. .. ..$ block.size        : num 1000
  .. .. .. .. ..$ min.cells.to.block: num 2638
  .. .. .. .. ..$ verbose           : logi TRUE
  .. ..$ RunPCA.RNA              :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
  .. .. .. ..@ name       : chr "RunPCA.RNA"
  .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2023-05-21 03:58:25"
  .. .. .. ..@ assay.used : chr "RNA"
  .. .. .. ..@ call.string: chr "RunPCA(seurat, npcs = 50)"
  .. .. .. ..@ params     :List of 10
  .. .. .. .. ..$ assay          : chr "RNA"
  .. .. .. .. ..$ npcs           : num 50
  .. .. .. .. ..$ rev.pca        : logi FALSE
  .. .. .. .. ..$ weight.by.var  : logi TRUE
  .. .. .. .. ..$ verbose        : logi TRUE
  .. .. .. .. ..$ ndims.print    : int [1:5] 1 2 3 4 5
  .. .. .. .. ..$ nfeatures.print: num 30
  .. .. .. .. ..$ reduction.name : chr "pca"
  .. .. .. .. ..$ reduction.key  : chr "PC_"
  .. .. .. .. ..$ seed.use       : num 42
  .. ..$ RunUMAP.RNA.pca         :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
  .. .. .. ..@ name       : chr "RunUMAP.RNA.pca"
  .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2023-05-21 03:58:35"
  .. .. .. ..@ assay.used : chr "RNA"
  .. .. .. ..@ call.string: chr "RunUMAP(seurat, dims = 1:20)"
  .. .. .. ..@ params     :List of 26
  .. .. .. .. ..$ dims                : int [1:20] 1 2 3 4 5 6 7 8 9 10 ...
  .. .. .. .. ..$ reduction           : chr "pca"
  .. .. .. .. ..$ assay               : chr "RNA"
  .. .. .. .. ..$ slot                : chr "data"
  .. .. .. .. ..$ umap.method         : chr "uwot"
  .. .. .. .. ..$ return.model        : logi FALSE
  .. .. .. .. ..$ n.neighbors         : int 30
  .. .. .. .. ..$ n.components        : int 2
  .. .. .. .. ..$ metric              : chr "cosine"
  .. .. .. .. ..$ learning.rate       : num 1
  .. .. .. .. ..$ min.dist            : num 0.3
  .. .. .. .. ..$ spread              : num 1
  .. .. .. .. ..$ set.op.mix.ratio    : num 1
  .. .. .. .. ..$ local.connectivity  : int 1
  .. .. .. .. ..$ repulsion.strength  : num 1
  .. .. .. .. ..$ negative.sample.rate: int 5
  .. .. .. .. ..$ uwot.sgd            : logi FALSE
  .. .. .. .. ..$ seed.use            : int 42
  .. .. .. .. ..$ angular.rp.forest   : logi FALSE
  .. .. .. .. ..$ densmap             : logi FALSE
  .. .. .. .. ..$ dens.lambda         : num 2
  .. .. .. .. ..$ dens.frac           : num 0.3
  .. .. .. .. ..$ dens.var.shift      : num 0.1
  .. .. .. .. ..$ verbose             : logi TRUE
  .. .. .. .. ..$ reduction.name      : chr "umap"
  .. .. .. .. ..$ reduction.key       : chr "UMAP_"
  .. ..$ FindNeighbors.RNA.pca   :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
  .. .. .. ..@ name       : chr "FindNeighbors.RNA.pca"
  .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2023-05-21 03:58:35"
  .. .. .. ..@ assay.used : chr "RNA"
  .. .. .. ..@ call.string: chr "FindNeighbors(seurat, dims = 1:20)"
  .. .. .. ..@ params     :List of 17
  .. .. .. .. ..$ reduction      : chr "pca"
  .. .. .. .. ..$ dims           : int [1:20] 1 2 3 4 5 6 7 8 9 10 ...
  .. .. .. .. ..$ assay          : chr "RNA"
  .. .. .. .. ..$ k.param        : num 20
  .. .. .. .. ..$ return.neighbor: logi FALSE
  .. .. .. .. ..$ compute.SNN    : logi TRUE
  .. .. .. .. ..$ prune.SNN      : num 0.0667
  .. .. .. .. ..$ nn.method      : chr "annoy"
  .. .. .. .. ..$ n.trees        : num 50
  .. .. .. .. ..$ annoy.metric   : chr "euclidean"
  .. .. .. .. ..$ nn.eps         : num 0
  .. .. .. .. ..$ verbose        : logi TRUE
  .. .. .. .. ..$ force.recalc   : logi FALSE
  .. .. .. .. ..$ do.plot        : logi FALSE
  .. .. .. .. ..$ graph.name     : chr [1:2] "RNA_nn" "RNA_snn"
  .. .. .. .. ..$ l2.norm        : logi FALSE
  .. .. .. .. ..$ cache.index    : logi FALSE
  .. ..$ FindClusters            :Formal class 'SeuratCommand' [package "SeuratObject"] with 5 slots
  .. .. .. ..@ name       : chr "FindClusters"
  .. .. .. ..@ time.stamp : POSIXct[1:1], format: "2023-05-21 03:58:36"
  .. .. .. ..@ assay.used : chr "RNA"
  .. .. .. ..@ call.string: chr "FindClusters(seurat, resolution = 1)"
  .. .. .. ..@ params     :List of 10
  .. .. .. .. ..$ graph.name      : chr "RNA_snn"
  .. .. .. .. ..$ modularity.fxn  : num 1
  .. .. .. .. ..$ resolution      : num 1
  .. .. .. .. ..$ method          : chr "matrix"
  .. .. .. .. ..$ algorithm       : num 1
  .. .. .. .. ..$ n.start         : num 10
  .. .. .. .. ..$ n.iter          : num 10
  .. .. .. .. ..$ random.seed     : num 0
  .. .. .. .. ..$ group.singletons: logi TRUE
  .. .. .. .. ..$ verbose         : logi TRUE
  ..@ tools       : list()

My objective is

Generate feature plot
Then label the cells based on only variable genes,right now I see all the genes I would like to only label based on top genes which I'm storing at top_gene but Im not sure how to incorporate it instead of every genes which I now see in my drop down menu
The after labelling cells i want to export the labelled cluster as csv file

What Now Im getting is html file which not what I want
Any help or suggestion is really helpful
UPDATE
The right function I have to use is this FetchData() but Im not sure how to pass the gene which i choose to label the cells and use the same gene as input for exporting cluster


text:
Evolutionary tree from SNP data
I would like to make an evolutionary tree using the SNP data from DArTseq. Are there any software or packages available for this purpose? I came across the software SNphylo, but its home page is not updated. Your suggestions would be appreciated.

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.

Welcome @Lama You would need to explain DArTseq and the data set size

At present I would simply use iqtree
iqtree -s file -m GTR+I+G -T 6 -B 1000

This will parallelise across 6 CPUs and results in 1000 bootstraps, whilst GTR is the mutation model, which models invariant sites under a 4 category gamma distribution. This is a robust model.
You can view the tree in FigTree.
The DNA file is in fasta format
> seq1
AGACGT
>seq2
AGAAAT

If there's an 'unknown' format bug let me know because I only work in  data pipeline so I forget if I've shifted the format.
The result will be in the treefile and is a professional quality tree, complete with bootstraps.
The I parameter here could omitted for this dataset because I assume all none of the sites are invariant.

@M_ Thank you for you answer. DArTseq is one of the cheap and easy but efficient genotyping-by-sequencing platforms which allow genome-wide marker discovery through restriction enzyme-mediated genome complexity reduction and sequencing of the restriction fragments (Edet et.al. 2018). I am attempting to construct a phylogenetic tree using the results of DArTseq rather than the FASTA sequences.

@Lama its a simple format shift. Maybe conversion of DArTseq to fasta is a good question to raise? At the end of the day the theory is universal if the output is nucleotide data. We'd need to see the new type first, but it's likely very easy. You can build trees from anything and everything - its the underlying theory that counts.


text:
How to add bootstrap values to the phylogenetic tree generated by OrthoFinder?
When we run the OrthoFinder analysis tool on a group of genomes to get the orthologues shared by them one of the output files include a folder named 'Species_Tree' that contains a text file named 'SpeciesTree_rooted.txt' . These are its contents (This is an example I have used for the question - 11 genomes of the Metabacillus genus) :

((Mlacus:0.244776,Mmangrovi:0.284609)1:0.0277466,(((Mdong:0.0284157,Midriensis:0.0250137)1:0.0553715,Mindicus:0.0904605)1:0.148105,(Mfastidiosus:0.207986,(Mniabensis:0.142138,((Mcrassostreae:0.0124653,Mlitoralis:0.0139468)1:0.116616,(Mhalosaccharovorans:0.0408989,Mschmidt:0.0422379)1:0.0872204)1:0.0402912)1:0.0849362)1:0.0750784)1:0.0277466);

I used iTOL to construct the tree using this file but I could only see the branch lengths and no way to add any bootstrap values if there had been any, to the tree as seen in most research papers - 99, 98, 100 etc. Any suggestions?

Bootstrapping is easy just use iqtree
iqtree -s file -m GTR+I+G -T 6 -B 1000

This will parallelise across 6 CPUs and results in 1000 bootstraps, whilst GTR is the mutation model, which models invariant sites under a 4 category gamma distribution. This is a robust model.
The file is aligned fasta format (its all data pipeline so I forget if I've shifted the format) - although this can be changed to any input format
The result will be in the treefile and is a professional quality tree, complete with bootstraps.
The following up analysis is a Bayesian tree and is more complicated.
Note You must resolve/delete gaps using iqtree because of the way the authors are calling gaps - this related to your question here - How to use Gblocks for trimming single-copy gene sequences? . Basically columns with gaps should be deleted.

Question, are you sure this is not a Bayesian analysis? It would make sense if it was a Bayesian analysis that is 100% robust (there are robustness scores there). This would make sense for this tree and this data. As far  I am aware iTOL is used for viewing and manipulating trees rather than creating them.

Thank you very much for the answer @M__ . Yes, iTOL is used for viewing trees and not creating them. The tree itself is created by Fasttree which is part of OrthoFinder, but the output is a text file. I did use this command and I got an error : File not found GTR .  And then I included -st DNA in the command but it still did not resolve the issue.  This is what I get -- WARNING: Number of command-line arguments differs from checkpoint
WARNING: Command-line differs from checkpoint!
Reading alignment file Metabacillus_Tree.fa ... Fasta format detected
...

Hi @KirtiKulanthaivel the best thing to do arrange a chat? If ok you will need your data available and access to a terminal and understand bash commands such as more and/or have access to a text editor such as BBEdit. GTR is central to the model, it will not work without it. Also could you install iqtree via conda? If that is not clear I will explain. If ok please specify availability in GMT, early afternoon is preferred. It is better to do this interactively than via question/answer. I will ask questions and ask you to post a small sample of your data and output errors. Let me know.

The bioinformatics chat room is here: https://chat.stackexchange.com/rooms/58858/bioinformatics. .. can alternatively do this via question/answer - the key thing is GTR must be resolved. Something is not right for this error to be thrown and I don't know what it is. The theory as to why is complicated but getting it working is getting it working.

Thanks @M__ how about 15:00 hours GMT?

@KirtiKulanthaivel my apologies, I am not around today, but tomorrow is fine. 15:00 hours cool.

Sure @M__ . See you!

Hi @M__ ! Just wanted you to know that the IQTREE code worked but I had to increase the bootstrap replicates from 100 to 1000 when the terminal read this message - #replicates must be >= 1000
The tree came out well, thank you very much. By the way, when running iqtree, there is a small table that shows the Gap/Ambiguity, Composition and p-value for the aligned SNP files in question (for 11 genomes in my case). It says 0.00 under Gap/Ambiguity and p-value and 'failed' under composition, although the tree is made. Why is that so?

Hi @KirtiKulanthaivel thanks for the feedback. The gap/ambiguity is good because you removed all "gaps" via trimal. The ambiguity, is the sequence is high quality - no ambiguities. The p-value of 'failed' I would need to look up. It sounds like it's attempting to estimate ambiguities/gaps according to a known distribution (e.g. Poison) but it can't because there aren't any. If it was a p-value <0.01 I suspect that would be a problem, 'failed' seems fine.


text:
Is there any python module to assign secondary structures which includes both Polyprolene region as well as individual strand as well?
I wanted to assign secondary structure to set of pdb files. However, DSSP and STRIDE cannot assign for single peptides if it is falling inside single long beta strand or polyprolene region as it does not have any hydrogen bond.DSSP and STRIDE work fine as long as the peptides form hydrogen bonds.
Can someone suggest me a way to assign secondary structure for peptides which also fall in the region of PPII and beta strand region ?


text:
How can I delete these lines in my fastq file?
Hi guys =) I'm sorry if this is a repeat but I haven't been able to find an answer using the search field or google.
I am trying to edit a fastq file that has a corrupt line. The read in question is VH00362:5:AAAV5KTHV:1:1210:64850:6869 where the read is just a massive repeat – not the expected 110bp. I am trying to remove two lines from this file and asking around  I've been told to use sed as in
gzip -c file.fastq.gz | sed "45500973,45500974d" > new.fastq.gz
but when I run wc -l  to check the number of lines the results are incoherent, as if it has deleted much more than the two lines I asked. Another issue that I'm seeing with this file is that every time I run gzip I get a invalid compressed data--format violated error, which tells me the file is corrupted but I'm doing my best to salvage whatever I can from the file as it is an irreplaceable sample.
Is there another way to simply delete these two lines? Any tool or command you could recommend? Any help is deeply appreciated.

Does gzip -k file.fastq.gz report an error? If so, I would first try to gunzip the file, or zcat file.fastq.gz > file.fastq and then operate on the decompressed one. I am not sure, but I suspect that a corrupted archive might not cat as expected and I would feel safer getting as much of the file decompressed as I can first. I am not even sure you will get anything, please try and report back.

This is not the thing I do these days but ...
gunzip -c file.fastq.gz | sed "45500973,45500974d" > new.fastq

You were using gzip instead of gunzip, and that the source of the error. Finally the out file isn't compressed, so you probably want:
gunzip -c file.fastq.gz | sed "45500973,45500974d" | gzip > new.fastq.gz

That should work (I tested it and yes it works).

A fastq file has four line-separated fields per sequence. If you suspect that a sequence is corrupted, then you would want to delete all four lines related to that sequence. If you just delete the first two lines (sequence identifier and raw sequence letters), the file will most probably throw an error down the road. In your case, the sequence would look something like this -
@VH00362:5:AAAV5KTHV:1:1210:64850:6869
GATTTGGGGTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

So, the code to delete these four lines from a fastq.gz file would be -
gunzip -c file.fastq.gz | sed "45500973,45500976d" | gzip > new.fastq.gz

This sed command specifies a range of lines and says that delete lines from line number 45500973 to 45500976.


text:
What is the current state-of-the-art for classifying the ancestry groups of VCF files?
I'm wondering which tool or program is:

Most updated; doesn't use a superseded algorithm
Runs very quickly (ideally, under 30s per VCF file)
Able to be automated in some way for a large number of VCFs
I need to at least have these groups: Africans, East Asians, Europeans, and South Asians. More specific groups are fine.
Able to be run locally, as I have a lot of data

I am not looking to visualize the data or graph a PCA. My goal is to get a percentage match, or any closeness measure, for each group for each VCF file (assuming one sample per file). All VCFs are generated from microarray data.
For example, Eurogenes k13 and MDLP World-22, especially with oracle, have sufficient results but are only available as a web app.
To clarify, I am asking if there exists any such local program that exists that is peer-reviewed or validated, fast, and automate-able.

Solution: use https://github.com/stevenliuyi/admix

Please add that as an answer.

I'm not sure that it's a good solution enough to post it because as far as I can tell, they're non-validated algorithms aimed at consumers, so it doesn't address the first bullet point.

Please mention / discuss your actual problem that you are trying to solve. General "what is the best" questions don't work well for Stack Exchange questions because there's no right answer that works in all cases.

If you're writing it in the comments as a solution to your question, it's good enough as the start of an answer. It can always be updated / improved at a later date.

'Runs very quickly (ideally, under 30s per VCF file)' - whats your sample size. One sample per vcf? #

Ah, good point. I really mean "per sample" instead of per VCF file. I was assuming one sample per vcf, yes.

Human genetics ain't my thing but I'd use a published machine learning or deep learning solution with an accuracy of >99%. This will be quick and accurate. It an obvious and easy supervised classification problem which means a large number of groups will already have done it.


text:
Off-target % for whole-exome sequencing panel
My samples have been sequenced with the Twist Exome 2.0 sequencing panel, and I want to assess how efficient the sequencing has been. By efficient, I mean examining metrics such as the uniformity (fold-80 score), coverage, on/off-target % etc. as documented in their webpage. I checked the % off-target as follows:

Count the # of reads in the bam file (samtools view -c "$dir/$base.bam)

Count the # of reads that intersect with the bed region (bedtools intersect -a "$dir/$base.bam" -b $bed -u | samtools view -c -)

Calculate the % off-target rate as (bc -l <<< "(1 - ($bed_reads / $total_reads)) * 100")

I have two questions:

Is there a cleaner/simpler way to calculate the off-target rate? I know I can calculate the on-target rate with bedtools coverage, but I am not sure of how to calculate the off-target rate with a single command line.
I obtained an off-target rate of approximately 20-25% for my samples. Is this within the "good" off-target % range?


text:
Getting VCF file that contain common SNPs from 6 VCF file using isec
I have 6 VCF files, where I would like to obtain the SNPs that are common (by position) in all the 6 files. I have tried this command
bcftools isec  -p trial2 -n=6 -c all A2_SNP.vcf.gz A_SNP.vcf.gz B_SNP.vcf.gz C_SNP.vcf.gzD_SNP.vcf.gz E_SNP.vcf.gz

But this produces 6 files, and through the README file, I think each of the VCF contains the unique SNPs for each of the 6 files, and I could only get sites.txt file that ocntain the positions that are common. Are there any way to get a VCF file that have the common SNPs ?
Thanks

You've used the right tool and gotten the right results. What you need to understand is that VCF contains variant loci as well as the genotypes for >=1 sample(s) at those loci. The intersect operation will extract common loci, but it does not know which genotypes you're interested in, which is why it generates one file per input, giving you the genotypes from each VCF for the common loci.
You could use the sites file to do a bcftools merge operation, limiting to just the common sites using the -R <regions-file> option so you'd get a VCF with all the samples across all 6 input VCFs but only for the variant loci.

Thanks a lot. For my understanding: Do you means that the output 6 files from isec tool I used is a VCF for each sample and contains the genotype (and other info) for only that common sites (that are also in site.txt) ? If so, do you assume that if I merge this using -R, then that merged VCF file contains all 6 samples. How can I get the genotype or alleles for the 6 individuals ?

Yes, each VCF output should contain information from its corresponding VCF input. You could either use sites file to merge all 6 original inputs or simply merge all 6 isec outputs get to the merged common-loci VCF.

So, merging the 6 VCF file will justput them side by side on  sample basis ? Could you advise on a code for that ?

Try the command, look at the output and come back if you have further questions.


text:
Should genotype imputation be ancestry specific?
I'm wondering if imputation, specifically Beagle, needs a reference panel that matches the sample's ancestry group. For example, Beagle documentation suggests the 1000 Genomes Project phase 3 reference panel, which contains all ancestry groups. Is this acceptable to use to impute the genome of, say, an East Asian, African, and European?

Basically, yes, you will get better results if you use a reference panel which contains samples with similar ancestry to your targets. Imputation works by matching haplotypes in your target samples to those in a reference panel. If the haplotypes in the target don’t exist in the reference panel, the variants won’t be imputed accurately.

We demonstrated that using TOPMed sequencing data as the imputation
reference panel improves genotype imputation quality in these
populations[Hispanic/Latino and African], which subsequently enhanced gene-mapping power for
complex traits. For rare variants with minor allele frequency (MAF) <
0.5%, we observed a 2.3- to 6.1-fold increase in the number of well-imputed variants, with 11-34% improvement in average imputation
quality, compared to the state-of-the-art 1000 Genomes Project Phase 3
and Haplotype Reference Consortium reference panels.

From this paper.

Is it still acceptable to use a general reference panel like 1000 Genomes Project Phase 3, or has everyone switched to ancestry-specific by now?

Most people use TopMed now using a remote server, or use their own custom ancestry specific one,  I don’t think many serious studies use 1000 genomes anymore

TopMed is similar to 1000 genomes in that they are multiethnic. I agree that TopMed is better for other reasons. However, the question of "does a multiethnic panel which includes some samples of ethnicity X improves accuracy for ethnicity X over a panel without ethnicity X?" is a different question than "if the panel consisted of only ethnicity X, would that improve accuracy over the multiethnic panel?"

Well yes as the abstract says, E.g. Latino ancestry individuals are imputed better using topmed because they it contains Latino samples. You can only get good imputation of non-common variants if you have individuals of matched ancestry in the reference panel.

Right. The study does a good job answering the first question, but doesn't address the second question specifically. (Which is where 1000 genomes might have an advantage if using an online imputation server as it allows for ancestry-specific imputation whereas TopMed doesn't.)

Why do you think topmed isn’t ancestry specific and 1000 genomes is? Actually what do you mean by ‘ancestry specific’ exactly anyway?

On online imputation servers, TopMed only has a multiethnic panel option. 1000 genomes allows you to select a single population, for example EAS, to use. Ancestry-specific here means data from a single ancestry group as opposed to data from multiple ancestry groups.

@BigMistake Jonathan Marchini did 1000g imputation. He observed that pooling all samples together achieved accuracy than imputing in individual populations independently. Don't split by population.

Thanks! This is what I was wondering. That definitely makes it easier. Seems like sample size outweighs the specificity of LD patterns for now. If you happen to have a link to where he tried that, I'd appreciate it.

@user172818 Any chance you have a paper or Blog I can refer to ? Our lab was pondering over the same question. Also based on what BigMistake said.

in my experience you should use a publicly cited reference genome like GRCh37 or GRCh38
https://www.sciencedirect.com/science/article/pii/S0888754317300058
https://gatk.broadinstitute.org/hc/en-us/articles/360035890951-Human-genome-reference-builds-GRCh38-or-hg38-b37-hg19

I use GRCh37 but I need to use a reference panel in addition to a reference genome

The question asks about a reference panel, not a reference genome.


text:
snakemake restricting resources for specific rule
I am trying to generate several snakemake rules. However, I would like to include some restrictions due to storage limitations
I first would like to copy the the bamfiles, 1 for each sample, to a staging area. However, while I have over 50 samples, I can only process three samples at a time.
Then I want to compare the bam file to chunked candidate files, 1 for each chromosome. This however, I would like to run on 23 cores per sample so that once each bam file is copied 69 cores are working.
finally, if all chunks for a given sample have been processed I want to clean up the staging area, copy the bamfile for the next sample and process the chunks for that sample.
My current implementation, see below, just continues copying bamfiles and crashes when I run out of strorage;
rule copy_bam:
    input:
        bamfile=os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam"),
        bamindex=os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam.bai")
    output:
        staging_bamfile="{staging_area}/{sample}.bam",
        staging_bamindex="{staging_area}/{sample}.bam.bai"
    wildcard_constraints:
        sample="[A-Za-z0-9_]+"
    resources:
        bam_copy=1
    shell:
        """
        cp {input.bamfile} {output.staging_bamfile}
        cp {input.bamindex} {output.staging_bamindex}
        """

rule process_chunks:
    input:
        chunkfile="{staging_area}/chunks/{chunk}_with_header.bcf.gz",
        bamfile=rules.copy_bam.output.staging_bamfile
    output:
        dynamic("{staging_area}/calls/{sample}/{chunk}.bcf")

    threads: 23
    shell:
        """
        process chunks
        """

rule cleanup:
    input:
        bamfile="{staging_area}/{sample}.bam",
        bamindex="{staging_area}/{sample}.bam.bai",
        processed=lambda wildcards: expand("{staging_area}/calls/{sample}/{chunk}.bcf", chunk=get_chunks()[wildcards.sample], staging_area=config['staging_area'], sample=wildcards.sample)
    output:
        touch("{staging_area}/{sample}/cleanup_done.txt")
    shell:
        """
        rm {input.bamfile}
        rm {input.bamindex}
        rm {input.processed}
        """

run with: snakemake --resources bam_copy=3 --cores 70

One thing to note about any approach here that uses Snakemake's resources is that they're defined per running job, versus with the "in use until something happens in a later job" idea here. Jame's Hawley's suggestion of temp() should help, just note that I don't think the resource tracking even matters in that case (the resource is only "in use" while copy_bam is actively running.)

One Snakemake feature you can use is temporary outputs.
If you make your outputs from the copy_bam rule temporary, they'll be cleared up after no other rule requires that file as an input.
Your Snakefile would look something like this:
rule copy_bam:
    input:
        bam = os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam"),
        idx = os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam.bai"),
    output:
        bam = temp("{staging_area}/{sample}.bam"),
        idx = temp("{staging_area}/{sample}.bam.bai"),
    wildcard_constraints:
        sample="[A-Za-z0-9_]+"
    resources:
        bam_copy=1
    shell:
        """
        cp {input.bam} {output.bam}
        cp {input.idx} {output.idx}
        """
# omitted for brevity

You can then specify how many jobs run concurrently with the --cores/--jobs options to keep the space usage small.


text:
Extracting information from VCF file for many specific positions in specific chromosomes?
I have an excel file that I created from VCF file for common SNPs across 6 samples. This excel have the chromosomes and the position of the SNPs only (see example table1)
Now I would like to obtain the other information (eg. allels, Genotype, depth, etc) from the VCF files of the 6 samples (i.e. the one that contains these positions).
I tried using AWK command like here for position 23432 on chr. 1 for the 6 file :
awk -F " " '$1=="1" && $2=="23432"' file1.vcf
awk -F " " '$1=="1" && $2=="23432"' file2.vcf
awk -F " " '$1=="1" && $2=="23432"' file3.vcf
awk -F " " '$1=="1" && $2=="23432"' file4.vcf
awk -F " " '$1=="1" && $2=="23432"' file5.vcf
awk -F " " '$1=="1" && $2=="23432"' file6.vcf

he issue is that these SNPs I have are thousands positions, so I need an automated way to do this
Could you advise on that ?
Thanks

Is this related to your previous question? Is this the sites file that you opened using Excel? Please read the bcftools manual, they have everything you need.

When working with VCF files, a tab-separated format, you want -F '\t'. And you never want -F " ", unless you truly want to break on spaces and only spaces. By default awk will split on any whitespace character (including tabs and spaces) so you don't need this.

Don’t use awk or excel, please just use bcftools and save yourself a lot of trouble

@user438383 OP has a separate thread going on with bcftools involved in solving this problem, looks like they chose to open a new thread for no apparent reason.


text:
StringTie discovers less novel isoforms when reference annotation is provided
I am analyzing PacBio IsoSeq data and I am using StringTie to assemble transcripts. I noticed the difference in the output files when I add a reference genome annotation (-G option in StringTie). In general I have less transcripts assembled when I give as imput reference genome annotation and I was wondering why. I know that the annotation of the genome is not the best (Borelian orangutang) and I would rather trust the data rather than the annotation provided.
I am attaching an example for one of the genes (IGV): the first track is reference annotation; the second is StringTie output when no reference annotation is provided; the third is StringTie output when reference annotation is provided (in the third case two novel assembled isoforms are missing)

Not an answer. For novel transcripts, you may consider IsoQuant. Its Nature Biotech paper argued that it does particularly better than StringTie on novel transcripts.

Not very familiar with what StringTie is doing, but in general reference annotations are helpful in providing a prior for what the tool expects.
I would suggest looking at the manual and using some of the provided options and outputs to probe the coverage of the various transcript species (I believe these are in the "attributes" section of the GTF output).
If the transcripts that are missing when you use -G have a coverage of just a couple of reads, they are likely to be less interesting. On the other hand, if they have relatively high coverage, then you have more reason to be skeptical of the existing annotation.
My guess is that the tool is more heavily weighting existing transcripts based on the annotation, and therefore excluding the extra transcripts as being below some threshold.
I'd recommend consulting the "recommended workflow" in the manual that supposedly allows for novel isoforms. Ensure that you are using this workflow and not the more streamlined workflow, which does not allow for novel isoforms.
See also the comment suggesting comparing to other IsoSeq-oriented tools.


text:
What do the numbers mean in these RNA-Seq gene/transcript TPM files?
From the link https://gtexportal.org/home/datasets, under V7, I'm trying to do R/Python analyses on the Gene TPM and Transcript TPM files. But in these files (and to open them I had to use Universal Viewer since the files are too large to view with an app like NotePad), I'm seeing a bunch of ID's for samples (i.e. GTEX-1117F-0226-SM-5GZZ7), followed by transcript ID's like ENSG00000223972.4, and then a bunch of numbers like 0.02865 (and they take up like 99% of the large files). Can someone help me decipher what the numbers mean, please? And are the numbers supposed to be assigned to a specific sample ID? (The amount of letters far exceed the amount of samples, btw). I tried opening these files as tables in R but I do not think R is categorizing the contents of the file correctly.
For context, I am planning to match males with females for sex comparison but in order to do that, I need to get R to categorize everything correctly. (I know that females have "F" where ####-#####-####-#x-##### where x is and males have "M").
Example subsets
For gene TPM:
#1.2 56202 11688
Name Description GTEX-1117F-0226-SM-5GZZ7 GTEX-111CU-1826-SM-5GZYN GTEX-111FC-0226-SM-5N9B8 GTEX-111VG-2326-SM-5N9BK GTEX-111YS-2426-SM-5GZZQ GTEX-1122O-2026-SM-5NQ91 GTEX-1128S-2126-SM-5H12U GTEX-113IC-0226-SM-5HL5C 7C-0005-SM-57WDL GTEX-ZVE2-0006-SM-51MRW GTEX-ZVZQ-0006-SM-51MR8 GTEX-ZXES-0005-SM-57WCB GTEX-ZXG5-0005-SM-57WCN
ENSG00000223972.4 DDX11L1 0.1082 0.1158 0.02104 0.02329 0 0.04641 0.03076 0.09358 0.121 0.02859 0 0.02373 0.08024 0.02231 0.0189 0.09309 0.07772 0.1596 0.1173 0.07969 0.07763 0.02062 0.07831 0.05159 0.1181 0.04187

For transcript TPM:
transcript_id gene_id GTEX-1117F-0226-SM-5GZZ7 GTEX-1117F-0426-SM-5EGHI GTEX-1117F-0526-SM-5EGHJ GTEX-1117F-0626-SM-5N9CS GTEX-1117F-0726-SM-5GIEN GTEX-1117F-1326-SM-5EGHH GTEX-1117F-2226-SM-5N9CH GTEX-1117F-2426-SM-5EGGH GTEX-111CU-0226-SM-5GZXC GTEX-111CU-0326-SM-5GZXO
ENST00000373020.4 ENSG00000000003.10 26.84 4.13 13.54 31.17 6.87 39.75 32.25 30.81 19.69 46.42 9.07 27.83 9.04 12.87 18.18 11.12 11.28 7.86 5.52 24.67 6.48 8.13 21.41 5.39 7.98 8.56 51.61 12.51 70.01 17.48 2.1 14.16 3.96 11.39 1.29 6.5 29.93 4.74 3.1 20.82 16.81 14.82 50.8 16.9 9.35 39.3 3.61 2.87 11.73

Is what I have there currently fine?

Thank you @Macromind101

The files have pretty self-explanatory headers.
The first file has EnsEMBL gene ID, HGNC symbol and then TPM for each sample. The second file has EnsEMBL gene ID, EnsEMBL transcript ID and TPM for each sample. Each header beginning with GTEX is a sample.
You should be using linux command line to look at this, not basic Windows utilities like Notepad. At least use Notepad++, basic Windows utilities are not equipped to handle large text files.

This is the first I've heard of linux command line. Is it a software application or a function in R? Did a quick google search and wasn't entirely sure.

Are you new to bioinformatics? Most of us use linux or macOS. You might be able to use WSL but without the command line, you're missing out on a lot of the power a computer has to offer in text processing. R and python are good but sometimes you need to do some quick operations and loading libraries/writing code to read a file can get tedious.


text:
Picard's CollectGcBiasMetrics gc-bias plot- help with interpretation
Setup
I ran picard's CollectGcBiasMetrics tool for my WES sample, with the intention of assessing the performance of the Exome capture panel in terms of the uniformity of read coverage at various %GC windows. I obtained the gc bias plot attached below.

Question
I am a bit puzzled by the normalised coverage values in the windows containing higher GC % and need help with interpretation - in particular, whether it indicates that something has gone awry with the sample, or if it is merely the way the normalised coverage is computed in picard. Why do the normalised coverage values seem to peak at around 50-55 GC%?
Digging into the gc_bias_metrics.txt file:

The #windows, #read starts and normalised coverage values are highlighted in the red, blue and yellow boxes respectively. Here, I notice that the # reads / # windows ratio is high for the tail end GC content, which probably inflates the normalised coverage value.
Is this something unexpected for a WES sample, i.e. should I be concerned?
Follow-up question
If the coverage values provided by picard isn't that useful, are there other useful tools for me to calculate the uniformity of read coverage across various GC%? I am looking at plots like Figure 3 in the twist biosciences data page.

Often coverage is higher in the relatively high GC range due to some artefactual biases from things like PCR, though the exome kits are getting better these days (and can be PCR-free): https://www.nature.com/articles/s41598-020-59026-y


text:
Is there any way to handle CONECT records in AutoDock Vina?
I modelled a macromolecule/receptor using Avogadro and exported it as receptor.pdb. Many entries in this file are CONECT records. If I try to run AutoDock Vina with it, this happens:

So, if I prepare this receptor again for docking, but, this time, removing the CONECT records (as I've seen many people do), my receptor file now exhibits broken bonds. Like this:

I've tried converting the .pdb file to .pdbqt with OpenBabel — removing all but the largest contiguous fragment and ignoring CONECT records as well —, but nothing works. I've also looked into it, and it looks like Vina doesn't have any parameters to handle CONECT records. Isn't there any way at all to do this?
Thank you in advance.

Hey, @M__! I'm sorry for not doing this sooner. I was a little short on time yesterday. What about now? Is there anything else I should add?

Hi, thats its better. If you copy and paste the ascii characters thats preferred to screen shots for the error output.