text:
single_to_multi_fast5 does not collect all the single files if the input folder contains mixed types of fast5 files
I have a dataset that contains thousands of mixed multiple and single fast5 files in a non-homogenous folder structure.
The reads and fast5 files are coming from different labs and some labs gave single fast5 files and some multi fast5 files, without clear mappings. My final goal is to convert all the fast5 files to multi fast5 files.
There are commands in ont_fast5_api which is an interface to HDF5 files of the Oxford Nanopore.
Using this API we can convert single fast5 files and multi fast5 files to each other.
this is the link to this API github page:
https://github.com/nanoporetech/ont_fast5_api
So I made a test mixed folder to check if I can use the API without missing reads as a result of conversion. The test folder structure is similar to my original folder structure:

mixed

Subfolder_A

single_fast5 x 4000

Subfolder_B

multiple_fast5 (contains 4000 single reads fast5)

single_fast5 x 2
multiple_fast5_1 (contains 4000 single reads fast5)
multiple_fast5_2 (contains 4000 single reads fast5)
multiple_fast5_3 (contains 100 single reads fast5)

Here is a screenshot of the structure with a few files from each folder:

and I created two other folders outside of the mixed folder as follows:

single

multi

Here my question is about single_to_multi_fast5 and multi_to_single_fast5 commands.
test_path=mixed
multi_path=multi
single_path=single
multi_to_single_fast5 -i $test_path/ -s $single_path/ --recursive

in the $single_path, I expect to see, all of the reads in multi files converted to single files (12100 reads in 12100 single fast5 files), but I see 8100 single files.

So I have 4000 missing reads.
Now I convert single read fast5 file to multi:
single_to_multi_fast5 -i $test_path/ -s $multi_path/ --filename_base $output_name --batch_size 1000 --recursive

in the $multi_path folder, I expect to see that generated multi files contain all my single-read files (4000 + 2), but I see 3006 in the mapping summary.

So again I missed 996 single reads in my multi files.
Since I made a small subsample of data, I am sure that the files are not overwritten.
Because my final goal is converting everything to multi I was going to do the following pipeline:
orig_path=original_folder
intermediate_path=intermediate
final_path=final_multi
# converts/collects all multi files in the original folder
multi_to_single_fast5 -i $orig_path/ -s $intermediate_path/ --recursive
single_to_multi_fast5 -i $intermediate_path/ -s $final_path/ --filename_base $output_name --batch_size 1000 --recursive

# converts/collects all single files in the original folder
single_to_multi_fast5 -i $orig_path/ -s $final_path/ --filename_base $output_name --batch_size 1000 --recursive

However because it seems I get missing reads during the process, I cannot use this pipeline.
At the same time, my original dataset is super huge, I cannot check if individual files are single or multi. It would take ages (I guess).
Is there a solution for solving this problem besides checking every file? Maybe my pipeline is not quite suitable for what I want to do.

Can you please show the first few filenames in each folder? e.g. by running the following command: find mixed -type 'd' | while read dname; do echo "** ${dname} **"; tree -L 1 ${dname} | head; done

Thank you for your reply! I included the first few filenames in each folder @gringer

The image seems like it's annotated incorrectly. A file with a channel and read number should be a single fast5 file, not a multi fast5 file.

Yes, the 2nd label (multi fast5 files with 4000 reads) and the 3rd label (single fast5 files) should be swapped. @gringer

Okay, that's good to know; thanks for the clarification.

It looks like you're creating too much additional work for yourself. What I'm trying to lead to with my comments is that the file names should help in working out if a file represents single reads or multiple reads, in which case the single reads can be shifted into a separate folder.
Based on the file names that I would expect, a search for files containing '_read_[0-9]' should be sufficient to fish out the single fast5 files:
find . -name '*_read_[0-9]*' > single_reads.txt
mkdir -p single_reads
cat single_reads.txt | while read fname; do mv -i ${fname} single_reads; done

[That single_reads.txt file could be split into multiple files if directory entries get too large.]
After which you should have a single folder that contains thousands of single-read fast5 files, and other folders that contain multi-fast5 files, removing the problem of needing to deal with mixed file types.

I am gonna try this, thank you! @gringer


text:
Krakenuniq mmap error
I built Krakenuniq on an HPC, and I've been trying to run it on my data using the following command:
krakenuniq --db DBDIR --threads 1 All_reads.fastq

This definitely isn't powerful enough to get any reasonable output, but this is just a test run. However, I got the following error:
classify: unable to mmap ./DBDIR/database.kdb: Invalid argument

I haven't been able to find anything about this error online. DBDIR is 1TB, which seemed huge, so I was wondering if there might have been a problem with the build, but I saw no issues in the log.

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.

Is DBDIR in your working directory? Did you perhaps mean to write --db ${DBDIR}?

That's not how it's done in the documentation, and in any case it causes its own error:

krakenuniq: unable to find --threads in $KRAKEN_DB_PATH (undefined) DBDIR appears to be in the working directory (ls lists it there).

I have run into the same issue with the default Kraken using a prebuilt database at the university HPC.

@schmiggle Did you ever find a solution to this? I can see the database file (for example, ./DBDIR/database.kdb) is there, but I get the same error.

I switched to using Kraken2 instead because of issues like this. Didn't ever find a solution.

source: Stack Exchange

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

sort: votes


text:
Parsing pre-2007 SMILES string
How would one parse the SMILES string
 BrC[2]:C[3]:C(:CH:CH:CH:@2):CH:CH:CH:CH:@3

I rely on tools like rdkit and OpenBabel to parse SMILES, but both tools aren't able to parse this string.
More specifically, this SMILES string comes from the supplementary information of a 2007 paper on solubility prediction. After doing some searching, it seems that this paper was received prior to the 2007 standardization of OpenSMILES (this might be why canonical parsers like rdkit have trouble with it).
After searching a bit more, I found an older PubChem sketcher that was able to convert the SMILES in question into a more recognizable format
 (BrC1=C(C(=C(C2=C1C(=C(C(=C2[H])[H])[H])[H])[H])[H])[H])

(parsing this was successful via rdkit, see molecule image below) but still need a way to convert similarly-formatted SMILES strings into the OpenSMILES standard/canonical format en mass.

Are the brackets dynamic? Thus one bracket pattern for one structure a different one for  a different structure or are they always as stated in the second expression? Hope that makes sense - it does make a difference to the answer.

I don't think they're dynamic but am not very familiar with the SMILES format specification.

"Welcome to Jurassic park" (theme music starts playing). Wow, that is proper ancient. And more similar to a SMARTS than a SMILES!
Your outputted version is kekulised and a machine-made mess —I don't think that's a good option. It ought to be something like [Br]c1c2c(ccc1)cccc2 if keeping your input order, but I think the canonical order is : [Br]c1ccc2ccccc2c1 (unchecked, simply typed out while counting carbons aloud). As you can see the standard isn't too dissimilar, so one could write a conversion schema on a push —although I am not sure if it's worth it.

@MatteoFerla Do you know of any software tools that can execute this conversion? If not, how would you go about converting the older format to a more standard one?

I have not encountered pre-2007 SMILES nor do I know of any tool to do so. Few sites would still exist 15 years later. C++ pre-11 generally compiles —pre-98 is the nightmare. Scripting modules may exists, but in Perl, not Python... So CPAN would be the place to search. But the question is whether this is worth it... after all your example molecule is bromonaphthalene.

I see, thanks. My problem is not parsing this particular molecule but rather a whole host of similarly-formatted ones (around 2k), so an automated parser would definitely help!

I assumed as much, I meant it is a rather small compound and there are hordes of catalogues. However, I had a look at your source and it has a name field. Why not just search the name with PubChem PUG API to get the SMILES?
PS. If you do find a Perl module in CPAN that does do the conversion or a similar converter do post the answer!

I couldn't find a conversion module, but the name field trick worked - thanks! Am posting more details as an answer.

I couldn't find a way to directly convert these SMILES strings to a more canonical format. However, thanks to @Matteo Ferla, I was able to get around the problem entirely by looking up the chemical names  and resolving the canonical SMILES programmatically.
More specifically, I used the CAS no column given in the supplementary information of the referenced paper to pull the SMILES via Chemical Identifier Resolver. See this answer for an explanation on how to do so.
EDIT: Some of the CAS nos couldn't be resolved, so I manually converted them using the PubChem sketcher tool mentioned in the question. There are quite a few (~160) of such molecules, so if anyone needs the full parsed dataset but doesn't want to spend the time manually parsing those molecules, let me know.

As it’s currently written, your answer is unclear. Please [edit] to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center.

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


text:
Select synonymous sites from a multiple sequence alignment
Could someone kindly recommend a tool or R package that can identify synonymous sites in a multiple sequence alignment?
I wish to select those taxa for tree reconstruction and other downstream analyses within subset of the alignment.
To respond to the questions below: yes, I am interested in identifying synonymous sites within a multiple sequence alignment (DNA) of ~70 virus samples with limited genetic diversity. My goal is to subset the multiple sequence alignment including synonymous sites only, because we believe that there are parallel non-synonymous mutations occurring that do not provide information on shared evolutionary history.
Thank you @M_ for the really helpful information below. In Mega, the S tool identifies singleton sites rather than synonymous sites in a multiple sequence alignment. I'm wondering if there is another tool or a PAML option to subset the synonymous sites?
Thank you in advance and best!

Do we also need to handle the translation, or is it already a protein alignment? Are we talking about two sequences or more? What kind of output do you need? Visual? Or text? Something else?

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.

I proposed Syn-SCAN (below). This was used for HIV and as you are looking at viruses seems a good fit.
There's a simpler way ... just make a invariable amino acid alignment, i.e. discard all amino acid variation. MEGA will do that, or use an algorithm. You then simply check the remaining nucleotide sequences for identity and discard the identity (lots of algorithms do this).
What you are left with is, is what you seek, PURE SYNONYMOUS mutations.

discard taxa with amino acid variation
retain taxa nucleotides with variation.

However, just discarding amino acid variation and retaining everything else will either by identity or synonymous. To be honest, its a virus so the amount of sequence identity between taxa is likely to be very little for protein loci.
Its really that simple. I honestly wouldn't over complicate it.
Sorry for the singleton vs. synonymous error, I haven't looked at MEGA in ages. Its a bit bizarre it discarded the only thing it was actually used for.

The short answer to the question is the 3rd codon position is the traditional approximation for synonymous mutation. A small proportion of 1st codon position sites are neutral and a small proportion of 3rd codon position are non-synonymous.
I agree with @terdon about and @Community regarding additional information sought, because synonymous mutations are heavily subject to compositional bias and that is a big subject. The key issue is the whether the taxa/species in the alignments are known for compositional bias and the genetic distance between different members of the alignment.
The traditional package used for synonymous mutations is/was MEGA, currently version 11. Albeit, recent 'upgrades' seem to remove the synonymous distance matrix calculation. HOWEVER, they may have transferred this capability onto its what appears to be upgraded GUI alignment viewer. What I think is if you:

load the alignment into the GUI;
click alignment (big button far left of menu);
you will now see the alignment;
in the middle of the menu bar are 4 buttons 'C', 'V', 'Pi' and 'S'
I think 'S' = synonymous
All synonymous sites in the alignment are now highlighted and there is likely a method of extracting this information.

C=constant, V=variable Pi=Nei approximation of diversity (population genetics)
If this is correct then those are the subset sought, please note however you will need to select/deselect the taxa manually and reassess the synonymous sites after each 'point and click' iteration (the S button will remain active I suspect between each selection/deselection). What you appear to want is a synonymous tree and this function I believe was removed from MEGA.
MEGA is rarely used except for teaching because its functionality is available elsewhere that permit a given algorithm to be integrated into a data pipeline, which MEGA would not permit (might have recently changed). Hence my recollection of MEGA is a little vague.

The package used for synonymous work is PAML and remains the undisputed standard without question in the context of formal selection analysis, which is heavily dependent on a estimate of synonymous mutations. What you refer to as 'synonymous mutation' isn't the same as the phylogenetic definition, i.e. your definition depends on the other members of the alignment, whereas this isn't an issue in the phylogenetic definition (because the tree circumvents this, i.e. delineates subsets) except if compositional bias occurs ... huge amounts of work is/has been dedicated to compositional bias.
Final point is the widely used phylogenetic library in R called ape will not perform a synonymous mutation calculation - I think, however phangorn library in R might do, but is more likely within the context of selection analysis. Thus to answer the question directly, I don't think R will do this calculation.

It transpires S= singleton (sorry) ....
Have a look at Syn-SCAN. The link is here https://hivdb.stanford.edu/pages/synscan.html . I have never used it and the download link is on the website. The good thing is its written in Perl so it would be easy to modify and check (if you can read/write Perl).
PAML: you can't detach the synonymous mutation calculation from the non-synonymous mutation calculation, without delving into the C-code (which I wouldn't recommend).

"Final point is whether R's ape": ape? Should that be API or something else entirely?

Thank you @M_ for the above! I edited my question in response. Thanks again!

Hi @terdon .. sorry ape is a popular phylogenetics library in R. Personally, I think its over-used. Its not an API ... I will clarify the text.

Thanks! As for "Personally, I think its over-used", that's how I feel about R :P


text:
Extract SNPs from multiple sequence alignment
I’m wondering if folks have recommendations for tools/scripts to extract SNP sites from a multiple sequence alignment of consensus bacterial genomes? Specifically, I am interested in a multiple sequence alignment containing variant (SNP) positions only, to be used as input for phylogeny-building programs including RAxML and BEAST.
I was previously using https://github.com/sanger-pathogens/snp-sites, but this is not being maintained and I’m running into errors that it seems many others are having (with regards to sequences being different lengths, although I have confirmed they are the same length).
Any recommendations would be great, thank you!

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.

Are you looking for identity, or are you variant calling, or finally perhaps looking a consensus sequence?

Yes, apologies for the lack of clarity. I have aligned consensus sequences, I'm hoping to extract SNP sites only into an alignment of SNP positions. Thank you!


text:
ATAC peak anatations
I am new in ATAC analysis.
I have a problem about how to transfer ATAC peak to gene symbol.
Here is the ATAC-seq: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM3271041
The GSM3271041_ATAC_sciCAR_A549_peak.txt.gz file contains five columns : (id, peak, chr, start, end). What I am doing is to use bedtools to transfer locations in chromosome to get gene symbols.
Here is an example from GSM3271041_ATAC_sciCAR_A549_peak.txt.gz
id,peak,chr,start,end
9,1-68591-69091,1,68591,69091
10,1-91105-91605,1,91105,91605

I use bedtools to get the gene symbol of id=9 , that is OR4F5.
I want to use these ''symbol-transfered'' ATAC-seq to integrate with RNA-seq which has gene symbols alreadly.
Am I right to do so ?

ATAC-seq data is fundamentally different from RNA-seq data because they measure different things.
RNA-seq data is typically mapped from raw sequencing reads to annotated gene transcripts.
This has a nice 1:1 correspondence because RNA is transcribed from genes.
ATAC-seq data, on the other hand, is used to measure accessible regions of the genome.
But accessible regions don't come from genes the same way that RNA comes from genes.
Accessible regions are scattered throughout the genome and most lie in intergenic regions where there is no annotation for a gene.
Because accessible regions don't have the same biological connection to genes that RNA transcripts do, you need to be more precise in how you want to pair up accessible regions and genes.
Are you looking for cis-regulatory elements (CREs) that regulate nearby genes?
Those will often be found within accessible regions, so you'd need additional data that connects an individual regulatory element to that gene.
Are you looking to create a catalogue of accessible regions within a certain cell type? Then you might want to look at ENCODE's candidate CREs database to see how your ATAC-seq peaks compares.
From there, you'll be able to find additional information about those peaks, such as whether other functional genomics data identifies these regions as active regulatory elements, which gene(s) these elements regulate, and more.
There are different options depending on what exactly you're trying to do, but there is no one way to transfer ATAC-seq peaks to gene IDs because there's no biological question there.

Thanks for your answer!
I think I want to map the peak in scATAC-seq to nearby genes.
I found this.
It's called peak annotation.

Annotating each with with the nearest gene(s) is a more specific question, that's great. Elsewhere on 10X's site you can see a description of their annotation process. Essentially, if you get a list of genes in a BED file, you can use bedtools closest with your peak list to find the nearest genes. The will usually work as a first pass for most applications, but like I said above, it will highly depend on your specific biological question.


text:
Questions about CreateGeneActivityMatrix function in Seurat
Seurat V3 has a function CreateGeneActivityMatrix, which can convert ATAC peaks to gene-level's matrix.
Here is the prototype:
CreateGeneActivityMatrix(
  peak.matrix,
  annotation.file,
  seq.levels = c(1:22, "X", "Y"),
  include.body = TRUE,
  upstream = 2000,
  downstream = 0,
  verbose = TRUE
)

Here are my questions:

I want to know what is the format of the first parameter peak.matrix ? For example, what are the row names of this matrix? (Is there a strict rule of the row names?)
The parameter seq.levels. Besides 1...22,X,Y , Can I add other chromosome ? like hs37d5


text:
Pulling Multiple Files via rsync with WDL
I'm trying to run a simple task to pull multiple files via rsync with WDL. The task should be quite simple but how to pull the multiple files subsequently is giving me issues. Below is a snippet of my task and how it's being called:
task MultiFiles {
  input {
    String wallclock
    String memory
    Int tmpspace_gb
    Int cpu
    Array[String] remote_files #files we want to pull
    Array[File] local_files
    String data_host
    Boolean? checksum_flag = false

  }
    #Optionality for WDL
    String checksum = if checksum_flag == true then "--checksum" else ""

  output {
    Array[File] output_files = "${local_files}"
  }

  command <<<

    set -e

    #rsync command
    rsync -avP
      --log-file=rsync.log \
      ~{checksum} \
      "~{data_host}":"~{remote_files}" \
      "~{local_files}"

  >>>
  runtime {
    cpu: 4
    memory: "2GB"
  }
}

Call
    call MultiFiles {
      input:
        #Task Resources
        remote_files=remote_files,
        data_host=data_host,
        local_files=local_files
}


text:
To find genes that don't change in RNA-seq, Deseq2 has altHypothesis="lessAbs". Is there a way to make limma do the same thing for proteomics?
I love that Deseq2 has altHypothesis="lessAbs" !!!!!!! I've used it a ton for RNA-seq. However, now I'm working with mass spect data using limma. Is there a way to make limma do the same thing?

Hi and welcome to the site! Could you please [edit] your question and clarify a bit? First, please make the title be a title: a short, concise description of what the question is about. Second, please give us a bit more context. I think limma is an R package, right? And what exactly would this option do? Identify genes that don't change? Between what and what? Are you comparing two sets of mass-spec data? Please add all extra info to your question, don't answer in the comments, just add it to the main question.

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.

There is no direct support for this by best knowledge. Similar questions have been asked previously in the Bioconductor support forum, for example:
https://support.bioconductor.org/p/64300/#64793
They discuss options around the TOST concept.
For edgeR:
https://support.bioconductor.org/p/66283/#66287
From these posts please go back and forth to find other posts in that forum to capture what the developers of limma/edgeR responded on the matter.
Depending on your exact purpose it might even be enough to go for genes with large AveExpr and a nominal (not adj) p-value close to one.


text:
Running cellranger on scRNASeq data with feature barcoding (10x + antibody capture)
I can't seem to find a clear answer to this question, so here it goes:
I have sequenced scRNASeq + scVDJSeq (TCR) data, which has been sequenced using feature barcoding from 10x genomics, via antibody capture.
There is a very useful cellranger multi command, which allows us to process both this data types in 'one go'.
In order to run it one needs to generate a csv file a priori, such that cellranger is able to identify which files are pertinent to which file type - see this link:https://support.10xgenomics.com/single-cell-vdj/software/pipelines/latest/using/multi
In the link above, with the example given under "Example multi config CSVs", it seems as though there are separate fastq files for the feature barcode, versus for the actual gene expression samples. Whereas, in my case, I only have the standard R1 and R2 fastq files, which is creating some confusion as how to generate this file, and run cellranger.
Now my question is, how does one run cellranger multi with feature barcoding should there only be the standard 2 fastq files? My understanding is that the feature barcode is usually in the R2 read? Is this correct? Would I need to Demultiplex? I am completely lost as to how to proceed. If you know, could you please provide a short example of what your csv files would look like, should you not have separate fastq files for the feature barcode files?
Thanks,

How did you create your fastq files? If you've run cellranger mkfastq with the correct sample sheet, you should have got separate fastq files for each library: RNA-seq, VDJ-seq and feature barcode.

The fastq files were provided by the sequencing company.. We have 4 fastq files (R1,R2,I1,I2) for the GEX & VDJ sequencing, all run via 10x genomics. i.e. I don't have any bcl files, just the fastq's + aligned BAM's

Before running cellranger multi you first need to demultiplex your fastq files to obtain separate files for each library. Unfortunately I've never demultiplexed fastq files so can't recommend any software.

Check your FASTQ headers - if they've been demultiplexed, the index barcodes can be found in the header. The barcodes are usually separated by a + in the comment part (second segment when the header is split by " " blank space) of the header. If you see a match to [ATGC]{10}[+][ATGC]{10} in the header, the first [ATGC]{10} is the I1 barcode and the second is the I2 barcode. You can use FFFFFFFFFF placeholder qual scores to generate "fake" I1 and I2 files and then use that with cellranger.

Thanks for your response. I'm not referring to the I1 & I2 files. I'm referring to feature barcoding technology using antibody capture (or potentially CRISPR). In the examples on the 10x website, they have separate fastq files for these, along with the gene expression fastq's

I see, I apologize for the misunderstanding - I do not know how to help you with feature barcoding.


text:
Looping over a list of cell models in Python
I have two data sets

Mutations in 10 genes
Gene expression in 18 cancer 3D models

For every (mutation, gene) and gene expression pair we run an ANOVA test
This is my data

> mut=read.delim("mut.tsv",row.names = 1)
> expr=as.data.frame(t(read.delim("expr.tsv",row.names = 1)))

dput(expr)

structure(list(GeneA_KO = c(0.0756383119639074, -0.162841965682723,
-0.110703100204517, -0.457154616239399, -0.505815193885352, 0.366072775033448,
0.288650458935083, -0.0287646498640309, -0.690732824311905, 0.157366122620995,
-1.29451546737684, 0.124167425460344, -0.168045313861751, 0.0186504943306465,
-0.495880827400706, -0.206668627339351, -0.148609490335447, 0.066844529089804
), GeneB_KO = c(-0.870157285953621, -0.35900872537626, -0.990406888128838,
-0.120656981579292, -1.41427631791049, -0.866421095836781, -0.816936112146099,
-0.894597236248507, -1.00084584207823, -0.850754488133456, -1.05314239761653,
-1.07580675428372, -1.24708324665688, -1.04098460460899, -0.232388675985086,
-1.48513703581865, -1.48373348870346, -1.76484058374354), class = "data.frame", row.names = c("ModelM",
"ModelO", "ModelS", "ModelA", "ModelI", "ModelC", "ModelT", "ModelH",
"ModelE", "ModelR", "ModelA2", "ModelP", "ModelE2", "ModelU",
"ModelT2", "ModelC2", "ModelS2", "ModelX"))

dput(mut)

structure(list(ModelM = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
0L), ModelO = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelS = c(0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelA = c(0L, 0L, 0L, 1L,
0L, 0L, 0L, 0L, 0L, 0L), ModelI = c(0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 0L), ModelC = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L
), ModelT = c(0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L), ModelH = c(1L,
1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L), ModelE = c(1L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 1L), ModelR = c(1L, 0L, 0L, 1L, 0L, 0L, 1L,
0L, 1L, 0L), ModelA2 = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
0L), ModelP = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L), ModelE2 = c(1L,
0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L), ModelU = c(1L, 0L, 1L, 1L,
0L, 1L, 0L, 1L, 1L, 1L), ModelT2 = c(0L, 0L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L), ModelI2 = c(1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L,
0L, 0L), ModelC2 = c(1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L),
    ModelS2 = c(0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L)), class = "data.frame", row.names = c("Gene1_mut",
"Gene2_mut", "Gene3_mut", "Gene4_mut", "Gene5_mut", "Gene6_mut",
"Gene7_mut", "Gene8_mut", "Gene9_mut", "Gene10_mut"))

To test this, for the first 3D cancer model as instance, I have done like below
           GeneA_KO   GeneB_KO
ModelM   0.07563831 -0.8701573
ModelO  -0.16284197 -0.3590087
ModelS  -0.11070310 -0.9904069
ModelA  -0.45715462 -0.1206570
ModelI  -0.50581519 -1.4142763
ModelC   0.36607278 -0.8664211
ModelT   0.28865046 -0.8169361
ModelH  -0.02876465 -0.8945972
ModelE  -0.69073282 -1.0008458
ModelR   0.15736612 -0.8507545
ModelA2 -1.29451547 -1.0531424
ModelP   0.12416743 -1.0758068
ModelE2 -0.16804531 -1.2470832
ModelU   0.01865049 -1.0409846
ModelT2 -0.49588083 -0.2323887
ModelC2 -0.20666863 -1.4851370
ModelS2 -0.14860949 -1.4837335
ModelX   0.06684453 -1.7648406

           ModelM ModelO ModelS ModelA ModelI ModelC ModelT ModelH ModelE
Gene1_mut       0      0      0      0      0      0      0      1      1
Gene2_mut       0      0      0      0      0      0      0      1      0
Gene3_mut       0      0      0      0      0      0      0      1      0
Gene4_mut       1      1      1      1      1      1      1      1      0
Gene5_mut       0      0      0      0      0      0      0      0      1
Gene6_mut       0      0      0      0      0      0      1      1      0
Gene7_mut       0      0      0      0      0      0      0      1      0
Gene8_mut       0      0      0      0      0      0      1      1      0
Gene9_mut       0      0      0      0      0      0      1      0      0
Gene10_mut      0      0      0      0      0      0      0      0      1
           ModelR ModelA2 ModelP ModelE2 ModelU ModelT2 ModelI2 ModelC2 ModelS2
Gene1_mut       1       0      0       1      1       0       1       1       0
Gene2_mut       0       0      0       0      0       0       0       1       1
Gene3_mut       0       0      0       0      1       0       0       0       0
Gene4_mut       1       1      1       1      1       1       1       0       1
Gene5_mut       0       0      0       0      0       0       0       0       1
Gene6_mut       0       0      0       1      1       0       0       1       1
Gene7_mut       1       0      0       0      0       0       0       0       1
Gene8_mut       0       0      0       0      1       0       0       1       0
Gene9_mut       1       0      0       0      1       0       0       0       0
Gene10_mut      0       0      1       0      1       0       0       0       1

> my_data=cbind(mut[,1],expr[,1])
> colnames(my_data)=c("mut","expr")
> res.aov <- aov(expr ~ mut, data = as.data.frame(my_data))
> # Summary of the analysis
> summary(res.aov)
            Df Sum Sq Mean Sq F value Pr(>F)
mut          1  0.145  0.1448   0.192  0.673
Residuals    8  6.022  0.7527

Now I do a for loop for all of the 18 cancer cell lines
> for(i in 1:ncol(g))
+ {
+   my_data=as.data.frame(cbind(m[,1],g[,1]))
+   #tidy will summarise and return neat format
+   print(broom::tidy(aov(my_data[,2] ~ my_data[,1], data = my_data)))
+ }
# A tibble: 2 × 6
  term            df sumsq meansq statistic p.value
  <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
1 my_data[, 1]     1 0.145  0.145     0.192   0.673
2 Residuals        8 6.02   0.753    NA      NA
# A tibble: 2 × 6
  term            df sumsq meansq statistic p.value
  <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
1 my_data[, 1]     1 0.145  0.145     0.192   0.673
2 Residuals        8 6.02   0.753    NA      NA

I done this in R but of course I got when doing in Python
Python
import os
import pandas as pd
from statsmodels.formula.api import ols
import statsmodels.api as sm
os.getcwd()
os.chdir("/Users//Downloads/test")
os.getcwd()
m = pd.read_csv("mut.tsv", sep="\t", index_col=0,delim_whitespace=True)
m=pd.DataFrame(m)
m
g = pd.read_csv("expr.tsv", sep="\t",index_col=0,delim_whitespace=True).T
g=pd.DataFrame(g)
g
lst = list(g.columns)

lst
dfs = []
for Model in lst:
   dfs.append(pd.concat([m.iloc[:, 0].reset_index(drop=True), g.iloc[0,:].reset_index(drop=True)],axis=1))
   dfs.append(df.dropna())
   dfs.append(df.columns==['mut','expr'])
   model = ols('expr ~ C(mut)', data=df).fit()
   aov_table = sm.stats.anova_lm(model, typ=2)
   aov_table
   final = pd.concat(dfs, ignore_index=True)

    for Model in lst:
      Input In [165]
        for Model in lst:
                         ^
    IndentationError: expected an indented block

When asking about errors, please remember to always include the actual error. We cannot help you debug errors you don't share.

thank you, I have added the error in my main post

@Angel, you still did not include the error that you get when you are using R. Anyway, why would you define i in your loop and not use it? And when it comes to the Python error, it is obvious!

Thank you @haci I get error in Python part only. The R code is OK

Thanks I have edited my post so I believe the R code will run smoothly. I have added the desire output from R code to the post. First, I have tried for one sample then using a for loop to doing one-way ANOVA over 18 samples. But, I do not know how to do this by Python

Here is one potential approach:
Example data:
expr1 <- structure(list(GeneA_KO = c(0.0756383119639074, -0.162841965682723, -0.110703100204517, -0.457154616239399, -0.505815193885352, 0.366072775033448, 0.288650458935083, -0.0287646498640309, -0.690732824311905, 0.157366122620995, -1.29451546737684, 0.124167425460344, -0.168045313861751, 0.0186504943306465, -0.495880827400706, -0.206668627339351, -0.148609490335447, 0.066844529089804), GeneB_KO = c(-0.870157285953621, -0.35900872537626, -0.990406888128838, -0.120656981579292, -1.41427631791049, -0.866421095836781, -0.816936112146099, -0.894597236248507, -1.00084584207823, -0.850754488133456, -1.05314239761653, -1.07580675428372, -1.24708324665688, -1.04098460460899, -0.232388675985086, -1.48513703581865, -1.48373348870346, -1.76484058374354)), class = "data.frame", row.names = c("ModelM", "ModelO", "ModelS", "ModelA", "ModelI", "ModelC", "ModelT", "ModelH", "ModelE", "ModelR", "ModelA2", "ModelP", "ModelE2", "ModelU", "ModelT2", "ModelC2", "ModelS2", "ModelX"))
mut <- structure(list(ModelM = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelO = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelS = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelA = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelI = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelC = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelT = c(0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L), ModelH = c(1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L), ModelE = c(1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L), ModelR = c(1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L), ModelA2 = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelP = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L), ModelE2 = c(1L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L), ModelU = c(1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L), ModelT2 = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelI2 = c(1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelC2 = c(1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L),     ModelS2 = c(0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L)), class = "data.frame", row.names = c("Gene1_mut", "Gene2_mut", "Gene3_mut", "Gene4_mut", "Gene5_mut", "Gene6_mut", "Gene7_mut", "Gene8_mut", "Gene9_mut", "Gene10_mut"))

R code:
my_data=cbind(mut[,1],expr1[1:10,1])
colnames(my_data)=c("mut","expr1")
res.aov <- aov(expr1 ~ mut, data = as.data.frame(my_data))
summary(res.aov)
#>             Df Sum Sq Mean Sq F value Pr(>F)
#> mut          1 0.1364  0.1364   1.114  0.322
#> Residuals    8 0.9789  0.1224

for (j in 1:nrow(mut)) {
  for(i in 1:ncol(expr1)) {
    my_data=as.data.frame(cbind(mut[,j],expr1[1:10,i]))
    print(broom::tidy(aov(my_data[,2] ~ my_data[,1], data = my_data)))
  }
}
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df  sumsq meansq statistic p.value
#>   <chr>        <dbl>  <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.0612 0.0612     0.465   0.515
#> 2 Residuals        8 1.05   0.132     NA      NA
#> # A tibble: 2 × 6
#>   term            df  sumsq meansq statistic p.value
#>   <chr>        <dbl>  <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.0637 0.0637     0.479   0.508
#> 2 Residuals        8 1.06   0.133     NA      NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.246  0.246      2.26   0.171
#> 2 Residuals        8 0.869  0.109     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.313  0.313      3.07   0.118
#> 2 Residuals        8 0.815  0.102     NA     NA
#> # A tibble: 2 × 6
#>   term            df   sumsq  meansq statistic p.value
#>   <chr>        <dbl>   <dbl>   <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.00108 0.00108   0.00777   0.932
#> 2 Residuals        8 1.11    0.139    NA        NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.220  0.220      1.94   0.201
#> 2 Residuals        8 0.907  0.113     NA     NA
#> # A tibble: 2 × 6
#>   term            df  sumsq meansq statistic p.value
#>   <chr>        <dbl>  <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.0529 0.0529     0.398   0.546
#> 2 Residuals        8 1.06   0.133     NA      NA
#> # A tibble: 2 × 6
#>   term            df  sumsq meansq statistic p.value
#>   <chr>        <dbl>  <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.0901 0.0901     0.695   0.429
#> 2 Residuals        8 1.04   0.130     NA      NA

# write.csv(mut, "mut.csv", quote = FALSE, row.names = FALSE)
# write.csv(expr1[1:10,], "expr1.csv", quote = FALSE, row.names = FALSE)

Created on 2023-09-19 with reprex v2.0.2

Python code:
import os
import pandas as pd
from statsmodels.formula.api import ols
import statsmodels.api as sm
m = pd.read_csv("mut.csv", sep=",")
m=pd.DataFrame(m)
m.info()
g = pd.read_csv("expr1.csv", sep=",")
g=pd.DataFrame(g)

lst = list(g.columns)
m_lst = list(m.columns)

dfs = pd.DataFrame({})
for Mut in m_lst:
    for Model in lst:
        df = pd.concat([m[[Mut]].reset_index(drop=True), g[[Model]].reset_index(drop=True)],axis=1)
        df = df.rename(columns={Mut: "mut", Model: "g"})
        model = ols('g ~ C(mut)', data=df).fit()
        aov_table = sm.stats.anova_lm(model, typ=2)
        aov_table[['Model']] = Model
        aov_table[['Mut']] = Mut
        dfs = dfs.append(aov_table)

with pd.option_context('display.max_rows', None, 'display.max_columns', None):
    print(dfs)
            sum_sq   df         F    PR(>F)     Model      Mut
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelM
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelM
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelM
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelM
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelO
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelO
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelO
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelO
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelS
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelS
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelS
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelS
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelA
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelA
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelA
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelA
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelI
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelI
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelI
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelI
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelC
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelC
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelC
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelC
C(mut)    0.061205  1.0  0.464539  0.514753  GeneA_KO   ModelT
Residual  1.054037  8.0       NaN       NaN  GeneA_KO   ModelT
C(mut)    0.063734  1.0  0.479234  0.508361  GeneB_KO   ModelT
Residual  1.063933  8.0       NaN       NaN  GeneB_KO   ModelT
C(mut)    0.245964  1.0  2.263619  0.170857  GeneA_KO   ModelH
Residual  0.869278  8.0       NaN       NaN  GeneA_KO   ModelH
C(mut)    0.312937  1.0  3.072788  0.117706  GeneB_KO   ModelH
Residual  0.814730  8.0       NaN       NaN  GeneB_KO   ModelH
C(mut)    0.001082  1.0  0.007771  0.931920  GeneA_KO   ModelE
Residual  1.114159  8.0       NaN       NaN  GeneA_KO   ModelE
C(mut)    0.220171  1.0  1.940909  0.201060  GeneB_KO   ModelE
Residual  0.907496  8.0       NaN       NaN  GeneB_KO   ModelE
C(mut)    0.052891  1.0  0.398296  0.545575  GeneA_KO   ModelR
Residual  1.062350  8.0       NaN       NaN  GeneA_KO   ModelR
C(mut)    0.090105  1.0  0.694741  0.428739  GeneB_KO   ModelR
Residual  1.037562  8.0       NaN       NaN  GeneB_KO   ModelR
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO  ModelA2
Residual  0.978877  8.0       NaN       NaN  GeneA_KO  ModelA2
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO  ModelA2
Residual  0.586718  8.0       NaN       NaN  GeneB_KO  ModelA2
C(mut)    0.004637  1.0  0.033399  0.859537  GeneA_KO   ModelP
Residual  1.110605  8.0       NaN       NaN  GeneA_KO   ModelP
C(mut)    0.276724  1.0  2.601573  0.145422  GeneB_KO   ModelP
Residual  0.850943  8.0       NaN       NaN  GeneB_KO   ModelP
C(mut)    0.044310  1.0  0.331000  0.580888  GeneA_KO  ModelE2
Residual  1.070932  8.0       NaN       NaN  GeneA_KO  ModelE2
C(mut)    0.170278  1.0  1.422852  0.267110  GeneB_KO  ModelE2
Residual  0.957389  8.0       NaN       NaN  GeneB_KO  ModelE2
C(mut)    0.001687  1.0  0.012120  0.915049  GeneA_KO   ModelU
Residual  1.113555  8.0       NaN       NaN  GeneA_KO   ModelU
C(mut)    0.008679  1.0  0.062048  0.809563  GeneB_KO   ModelU
Residual  1.118988  8.0       NaN       NaN  GeneB_KO   ModelU
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO  ModelT2
Residual  0.978877  8.0       NaN       NaN  GeneA_KO  ModelT2
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO  ModelT2
Residual  0.586718  8.0       NaN       NaN  GeneB_KO  ModelT2
C(mut)    0.017610  1.0  0.128353  0.729416  GeneA_KO  ModelI2
Residual  1.097631  8.0       NaN       NaN  GeneA_KO  ModelI2
C(mut)    0.260821  1.0  2.407079  0.159387  GeneB_KO  ModelI2
Residual  0.866846  8.0       NaN       NaN  GeneB_KO  ModelI2
C(mut)    0.191207  1.0  1.655405  0.234206  GeneA_KO  ModelC2
Residual  0.924035  8.0       NaN       NaN  GeneA_KO  ModelC2
C(mut)    0.033474  1.0  0.244741  0.634101  GeneB_KO  ModelC2
Residual  1.094193  8.0       NaN       NaN  GeneB_KO  ModelC2
C(mut)    0.044621  1.0  0.333425  0.579533  GeneA_KO  ModelS2
Residual  1.070620  8.0       NaN       NaN  GeneA_KO  ModelS2
C(mut)    0.096956  1.0  0.752533  0.410940  GeneB_KO  ModelS2
Residual  1.030711  8.0       NaN       NaN  GeneB_KO  ModelS2
```

Sorry I was running your python code but I get this error for Mut in m_lst:
Input In [81]
for Mut in m_lst:
^
IndentationError: expected an indented block

The IndentationError means you haven't got the correct whitespace - you need 4 spaces in front of each indented line @Angel and no whitespace in front of 'non-indented' lines (basically, copy it exactly as I have it in my answer and it will work)

Thanks a lot that's weird I just copied and pasted that but gives error

You're using spyder, so perhaps try going to "Source --> Fix Indentation" and see if your code changes and if the error persists

Source fixed Indentation error now I get this Traceback (most recent call last):
File "/Users/angel/Downloads/test/test.py", line 20, in
model = ols('g ~ C(mut)', data=df).fit()
File "/Users/angel/Downloads/docker/.venv/lib/python3.9/site-packages/statsmodels/base/model.py", line 209, in from_formula
raise ValueError('endog has evaluated to an array with multiple '
ValueError: endog has evaluated to an array with multiple columns that has shape (10, 18). This occurs when the variable converted to endog is non-numeric (e.g., bool or str).

That error is harder to understand, sorry @Angel; I think you should post another question with a reproducible example showing how you got that error so we can figure out what the problem is.

Before using this-->  with pd.option_context('display.max_rows', None, 'display.max_columns', None):,  print(dfs) could display the 1st 10 rows,  tried, worked!

Moreover, after using this--> python -V,  my z shell on my macbook air is using Python 3.9.16 .

Apparently the pd.option_context('display.max_rows', None, 'display.max_columns', None) approach doesn't work with jupyter notebooks; see https://stackoverflow.com/questions/16424493/pandas-setting-no-of-max-rows for more details @envs_h_gang_5


text:
Converting a data frame to multiple lists
I have data frame each row is a pathway and it involved genes like

like
> head(h)
    Release.of.Cancer.Cell.Antigens.1  ARID1A    ATM  BIRC3   BLM BNIP3   BRCA1   BRCA2 BRD3  BRD4 BRIP1 CCNA1 CCND1
V1        Cancer Antigen Presentation    ARG2   ATF3    B2M BATF3   BLK   CCL22 CCL3/L1 CCL4  CCL5  CCR5  CD14  CD19
V2      T-cell Priming and Activation ADORA2A  BATF3   BTLA  CCL2 CCL22 CCL3/L1    CCL4 CCL5  CCR2  CCR5  CD1C   CD2
V3 Immune Cell Localization to Tumors     A2M ADGRE1 ANGPT2  APOE  ATF3     AXL   BATF3  BLK CCL13 CCL14 CCL18 CCL19

I want to convert that like
gene.sets <- list(Release.of.Cancer.Cell.Antigens = c("ARID1A", "ATM",  "BIRC3",    "BLM"),
            Cancer_Antigen_Presentation = c("ARG2","ATF3","B2M"))

I have not done anything so far as I am not certain what exactly I should google for.
By code the user kindly provided I get this
gene.sets = lapply(1:nrow(h), function(x) unname(unlist(h[x,-1])))
names(gene.sets) = h[,1]

[91] "IGF2R"     "IHH"       "IL11"      "IL11RA"    "IL12RB2"   "IL15"      "IL17A"     "IL1A"      "IL1B"
[100] "IL22RA1"   "IL2RA"     "IL2RB"     "IL6"       "IL7R"      "IRF1"      "IRF2"      "IRF3"      "IRF4"
[109] "IRF7"      "IRF9"      "ISG15"     "ITGAE"     "JAK1"      "JAK2"      "JAK3"      "KIT"       "LAG3"
[118] "LAIR1"     "LIF"       "LTB"       "MX1"       "NCAM1"     "NKG7"      "OAS1"      "OAS2"      "OAS3"
[127] "PDCD1"     "PDCD1LG2"  "PF4"       "PIK3CA"    "PIK3CG"    "PIK3R5"    "PRF1"      "PRLR"      "PROM1"
[136] "PSMB10"    "PSMB5"     "PSMB9"     "PTEN"      "REN"       "RORC"      "RSAD2"     "SELL"      "SELP"
[145] "SIRPA"     "SOX11"     "SOX2"      "SPRY4"     "STAT1"     "STAT2"     "STAT4"     "TBX21"     "TIGIT"
[154] "TNF"       "TNFAIP3"   "TNFRSF10B" "TNFRSF10C" "TNFRSF11A" "TNFRSF11B" "TNFRSF14"  "TNFRSF17"  "TNFRSF18"
[163] "TNFRSF1A"  "TNFRSF1B"  "TNFRSF4"   "TNFRSF8"   "TNFRSF9"   "TNFSF10"   "TNFSF12"   "TNFSF13"   "TNFSF13B"
[172] "TNFSF18"   "TNFSF4"    "TNFSF8"    "TPSAB1/B2" "TRAF1"     "TSLP"      "VEGFA"     "WNT7B"     ""
[208] ""          ""          ""          ""          ""          ""          ""          ""          ""
[217] ""          ""          ""          ""          ""          ""          ""          ""          ""
[226] ""          ""          ""          ""          ""          ""          ""          ""          ""
[235] ""          ""          ""          ""          ""          ""          ""          ""          ""
[244] ""          ""          ""          ""          ""          ""          ""          ""          ""
[253] ""          ""          ""          ""          ""          ""          ""          ""          ""

x=lapply(gene.sets, function(z){ z[!is.na(z) & z != ""]})

I don't know whether this will actually work since you didn't provide a reproducible example...
gene.sets = lapply(1:nrow(h), function(x) unname(unlist(h[x,-1])))
names(gene.sets) = h[,1]

Thanks worked. As different pathways have different number of genes involved to I got a lot of empty thing. The pathway with the most number of genes has 292 genes so the code assumes each pathway should have 292 genes that is why for the rest puts "". I have added the result of your code in my main post

Thanks again I have added solution to remove ""


text:
Help Understanding the Greedy Motif Search Algorithm from Textbook
I am trying to go through this text book on Bioinformatics: https://www.bioinformaticsalgorithms.org. It has coding challenges and problems which I hadn't had too much trouble with until I got to the greedy motif search problem.
It is described in the book as such:

Our proposed greedy motif search algorithm, GreedyMotifSearch, starts by forming a motif matrix from arbitrarily selected k-mers in each string from Dna (which in our specific implementation is the first k-mer in each string). It then attempts to improve this initial motif matrix by trying each of the k-mers in Dna1 as the first motif. For a given choice of k-mer Motif1 in Dna1, it builds a profile matrix Profile for this lone k-mer, and sets Motif2 equal to the Profile-most probable k-mer in Dna2. It then iterates by updating Profile as the profile matrix formed from Motif1 and Motif2, and sets Motif3 equal to the Profile-most probable k-mer in Dna3. In general, after finding i − 1 k-mers Motifs in the first i − 1 strings of Dna, GreedyMotifSearch constructs Profile(Motifs) and selects the Profile-most probable k-mer from Dnai based on this profile matrix. After obtaining a k-mer from each string to obtain a collection Motifs, GreedyMotifSearch tests to see whether Motifs outscores the current best scoring collection of motifs and then moves Motif1 one symbol over in Dna1, beginning the entire process of generating Motifs again.

Pseudocode is also given:
GreedyMotifSearch(Dna, k, t)
    BestMotifs ← motif matrix formed by first k-mers in each string from Dna
    for each k-mer Motif in the first string from Dna
        Motif1 ← Motif
        for i = 2 to t
            form Profile from motifs Motif1, …, Motifi - 1
            Motifi ← Profile-most probable k-mer in the i-th string in Dna
        Motifs ← (Motif1, …, Motift)
        if Score(Motifs) < Score(BestMotifs)
            BestMotifs ← Motifs
    return BestMotifs

With this note:

Input: Integers k and t, followed by a collection of strings Dna.
Output: A collection of strings BestMotifs resulting from applying GreedyMotifSearch(Dna, k, t). If at any step you find more than one Profile-most probable k-mer in a given string, use the one occurring first.

I've exhausted other resources trying to understand this at the moment. The Output description basically says "the output is the strings which are output by this function" which is less than helpful. I don't understand what the 't' parameter is for, other than what appears to be an arbitrary iteration stopping point.
"For a given choice of k-mer Motif1 in Dna1, it builds a profile matrix Profile for this lone k-mer, and sets Motif2 equal to the Profile-most probable k-mer in Dna2." is the sentence where the text completely loses me. This is equivalent section in the pseudocode.
for each k-mer Motif in the first string from Dna
        Motif1 ← Motif

What is a k-mer motif? As far as I know, so far in the text, we have seen k-mers, we have seen them used to find motifs, we have seen profile-most probable k-mers, but at no point has anything been described as a k-mer motif. If someone has already gone through this book who could break down what it means in more detail it would be much appreciated.
EDIT: Just a note if anyone wants to tackle this problem and isn't already going through the textbook, the chapter in the book that I linked is free to access (as are the other 4 of the first 5 chapters in the book) and the section with this particular coding challenge is chapter 2.5: Greedy Motif Search. [ Lesson 2.5 ] (https://www.bioinformaticsalgorithms.org/bioinformatics-chapter-2)

I figured this out after rereading the chapter a few times. I am going to try my best to describe the solution here for anyone else who runs into this issue, since I don't think this was meant to be as difficult for the reader as it is, or at least as it was for me, especially since the function is described as "simple" in this video by the authors: Author's video on this section of the book. The part on GreedyMotifSearch is right at the end and lasts roughly a minute. Firstly, to clear up the input and output information:

Input: Integers k and t, followed by a collection of strings Dna. Output: A collection of strings BestMotifs resulting from applying GreedyMotifSearch(Dna, k, t). If at any step you find more than one Profile-most probable k-mer in a given string, use the one occurring first.

In defense of the writers of this, it had been a while since I had completed chapter 1 of this textbook, so maybe I just forgot the formatting that they had been using already. Regardless, I think the output explanation is still mostly void of helpful context.
Input:

dna: a list of dna strings to find motifs in
k: the character size of the motifs to be found in each string in dna
t: the number of strings that will be passed into the function. Maybe this is an easy way to do things in other coding languages, but I am using python. For me, it was more readable to remove this parameter and just reference the length of the dna parameter when I needed it using len(dna)

Output:

This is a list of motifs, one for each string passed in. The idea is to minimize the score of a list of kmers from each string in dna, and return the kmers that resulted in that minimum score. For me, Gringer's comment above helped a lot to clarify this.

Here is my version of the pseudocode:
GreedyMotifSeach(dna, k)
    best_motifs <- arbitrarily set to the first kmer in each string in dna
    for each 'kmer' in the first string in dna
        potential_better_motifs <- one item list containing 'kmer'
        potential_better_motifs_profile <- generate_profile_matrix(potentially_better_motifs)
        for string in len(dna) #  except for the first string. This is where the original pseudocode uses the t parameter.
            potential_better_motifs <- append the result of find_most_probable_kmer(dna[i], k, potentially_better_motifs_profile)
            potential_better_motifs_profile <- generate_profile_matrix(potential_better_motifs)
    if generate_score(potential_better_motifs) < generate_score(best_motifs)
         best_motifs <- potential_better_motifs
return best_motifs

So the part that I got lost at from my initial post:

"For a given choice of k-mer Motif1 in Dna1, it builds a profile matrix Profile for this lone k-mer, and sets Motif2 equal to the Profile-most probable k-mer in Dna2." is the sentence where the text completely loses me. This is equivalent section in the pseudocode.

Described the setting up of the potential_better_motifs and potential_better_motifs_profile variables which are changed throughout the internal loop, which goes through all of the strings in dna except for the first one.
For each kmer in the first string, a profile is built, so we always start with a profile matrix of 1s and 0s, since there is only one kmer that the profile is calculated from. Based on that profile a kmer is pulled from the second string. Then a profile is built from the those two kmers and from that a kmer is pulled from the third string. Rinse and repeat this process until you have a set of kmers, or motifs, from all strings in dna. This is greedy in the sense that the algorithm is only looking one move ahead, with only one kmer from string 1 being used to form a profile that finds the best kmer in string 2, only 2 strings being used to form a profile that finds the best kmer in string 3, etc.
With a new list of potentially_better_motifs, you need a way to compare your motifs to the best_motifs. To do this, generate a consensus string from the potentially_better_motif_profile, then count how many times that each string in your potentially_better_motifs differ from that consensus string. That count is your score for your potentially_better_motifs. Repeat this scoring process for your best_motifs, then compare the two scores. If the new motif set has a score that is lower (since this is a minimization problem, lower score == better consensus match to the profile), then you have a new set of best_motifs.
Repeat this whole process for each kmer in the first string in dna, then the function ends.
There are a lot of helper functions that I don't have pseudocode for. I think the book does a good job at describing how to implement these, and if I were to go through them all I would be coming close to rewriting the whole chapter in this post, but these are the helper functions that I used and their parameters.

find_most_probable_kmer(string, k, profile_matrix)
generate_count_matrix(kmers)
generate_profile_matrix(kmers)
generate_consensus_string(profile_matrix)
generate_score(kmers)

I chose to have generate_count_matrix be a helper for generate_profile_matrix. On top of using generate_profile_matrix within the GreedyMotifSearch, will likely be useful for you in generate_score. The parameters I used for most of these are just how I went about implementing this, so there is probably a more efficient method, but the parameters for find_most_probable_kmer are directly from the book, since that was the previous coding challenge. It is definitely worth noting the exact parameters for it, since those went a long way in helping me understand the core for loop in the GreedyMotifSearch function.

I have been similarly stuck on the Bioinformatics textbook motif problems in Rosalind, mostly because the documentation included with the problems is sparse, useless, or nonexistent.
The mouseover for "motif" in the SIMS problem calls a motif "a nucleotide or amino acid pattern of biological significance". In other words, a sequence/pattern that is being searched for inside another sequence.

What is a k-mer motif?

The pseudo-code doesn't have any punctuation to help understand it. In the code, Motif is the name of the variable used in the loop, and k-mer describes the elements of the sequence array.
In other words, "For each k-mer, Motif, in the first string from DNA..."
Consider the example input presented in the Rosalind GreedyMotifSearch problem (which should look familiar):
3 5 # kmer size of 3; 5 strings within the set of strings called "DNA"
GGCGTTCAGGCA
AAGAATCAGTCA
CAAGGAGTTCGC
CACGTCAATCAC
CAATAATATTCG

So in the first string from DNA, GGCGTTCAGGCA, the set of k-mers that represents possible motifs is [GGC, GCG, CGT, GTT, TTC, TCA, CAG, AGG, GGC, GCA]. These are processed as the variable Motif in the loops of the GreedyMotifSearch pseudo-code.


text:
BAM files with no RNAME and POS, how to map contents to SNPs?
I have a set of 4 .bam files containing the exome of an individual, around 400 MB each. I used samtools to generate a 2.4 GB .sam file out of one of the .bam files, and I found it contains lines with 101 base sequences but seemingly no RNAME or POS information on them. The first columns of the first lines in the read like this (I won't include here the SEQ, QUAL and optional fields):
210713_A00776_0149_AHCN5JDSX2:1:1101:10004:18818    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10004:18818    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10004:32033    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10004:32033    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:12712    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:12712    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:21261    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:21261    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:34882    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:34882    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:07764    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:07764    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10031:11569    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10031:11569    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10031:26537    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10031:26537    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10041:21010    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10041:21010    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10041:27085    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10041:27085    141 *   0   0   *   *   0   0

Most fields are nullified. As far as I can understand, subsequent lines with the 77 and 141 flags form corresponding pairs of bases and the QNAME field first element is an individual identifier, as it won't change through the head and tail of the file. Same goes for the second element (1). The last three elements change and might hold the information I need, but I'm not sure how to read them.
What I want to do is to map bases to SNPs (or a nucleotide identifier like NC_000006.11).
I have no real experience with bioinformatics besides enriching a file from MyHeritage (which contains around 609,000 base pairs with a corresponding rSNP number and chromossome) and crossed it with a reference human genome database I downloaded from NBCI. As it's a personal project, I somewhat validated the final result by correctly spotting a few known mutations/polymorphisms for that individual tested through usual lab tests, like the one output by the code I've wrote:
 RS1799945 (Chromosome 6 ): CG (NC_000006.11:g.26091179C>G, snv), https://www.ncbi.nlm.nih.gov/snp/rs1799945
          Reference Allele: C
                 Frequency:  9.6% (0.09631744626441) (~ 1:10)
              Loci (Genes): HFE  (homeostatic iron regulator)
                          : LOC108783645  (HFE antisense RNA)
        Detected Mutations: c.187C>G
             Variation URL: https://www.ncbi.nlm.nih.gov/clinvar/variation/10/
                  Diseases: not provided (conflicting-interpretations-of-pathogenicity), None (other), Hereditary hemochromatosis (pathogenic), Cystic fibrosis (risk-factor), Alzheimer disease (pathogenic), Familial porphyria cutanea tarda (None), Hemochromatosis type 1 (None), Microvascular complications of diabetes 7 (None), Transferrin serum leve, Microvascular complications of diabetes 7 (risk-factor), Hemochromatosis type 1 (pathogenic), Hemochromatosis type 2 (pathogenic)
                    PubMed: https://www.ncbi.nlm.nih.gov/pubmed/8696333, https://www.ncbi.nlm.nih.gov/pubmed/8896550, https://www.ncbi.nlm.nih.gov/pubmed/8943161, https://www.ncbi.nlm.nih.gov/pubmed/9106528, https://www.ncbi.nlm.nih.gov/pubmed/9326341, https://www.ncbi.nlm.nih.gov/pubmed/9462220, https://www.ncbi.nlm.nih.gov/pubmed/10673304, https://www.ncbi.nlm.nih.gov/pubmed/11358905, https://www.ncbi.nlm.nih.gov/pubmed/11423500, https://www.ncbi.nlm.nih.gov/pubmed/11812557, https://www.ncbi.nlm.nih.gov/pubmed/11874997, https://www.ncbi.nlm.nih.gov/pubmed/11904676, https://www.ncbi.nlm.nih.gov/pubmed/12429850, https://www.ncbi.nlm.nih.gov/pubmed/14673107, https://www.ncbi.nlm.nih.gov/pubmed/15347835, https://www.ncbi.nlm.nih.gov/pubmed/15858186, https://www.ncbi.nlm.nih.gov/pubmed/18194558, https://www.ncbi.nlm.nih.gov/pubmed/18199861, https://www.ncbi.nlm.nih.gov/pubmed/18566337, https://www.ncbi.nlm.nih.gov/pubmed/18603647, https://www.ncbi.nlm.nih.gov/pubmed/18795173, https://www.ncbi.nlm.nih.gov/pubmed/19084217, https://www.ncbi.nlm.nih.gov/pubmed/19159930, https://www.ncbi.nlm.nih.gov/pubmed/19165391, https://www.ncbi.nlm.nih.gov/pubmed/19237423, https://www.ncbi.nlm.nih.gov/pubmed/19401444, https://www.ncbi.nlm.nih.gov/pubmed/19554541, https://www.ncbi.nlm.nih.gov/pubmed/19862010, https://www.ncbi.nlm.nih.gov/pubmed/19879291, https://www.ncbi.nlm.nih.gov/pubmed/19884647, https://www.ncbi.nlm.nih.gov/pubmed/20110814, https://www.ncbi.nlm.nih.gov/pubmed/20301613, https://www.ncbi.nlm.nih.gov/pubmed/20435227, https://www.ncbi.nlm.nih.gov/pubmed/20556870, https://www.ncbi.nlm.nih.gov/pubmed/20659343, https://www.ncbi.nlm.nih.gov/pubmed/20846924, https://www.ncbi.nlm.nih.gov/pubmed/20876667, https://www.ncbi.nlm.nih.gov/pubmed/21138591, https://www.ncbi.nlm.nih.gov/pubmed/21149283, https://www.ncbi.nlm.nih.gov/pubmed/21208937, https://www.ncbi.nlm.nih.gov/pubmed/21540342, https://www.ncbi.nlm.nih.gov/pubmed/21553046, https://www.ncbi.nlm.nih.gov/pubmed/21679129, https://www.ncbi.nlm.nih.gov/pubmed/21909115, https://www.ncbi.nlm.nih.gov/pubmed/21978626, https://www.ncbi.nlm.nih.gov/pubmed/22074419, https://www.ncbi.nlm.nih.gov/pubmed/22100073, https://www.ncbi.nlm.nih.gov/pubmed/22125427, https://www.ncbi.nlm.nih.gov/pubmed/22232660, https://www.ncbi.nlm.nih.gov/pubmed/22408404, https://www.ncbi.nlm.nih.gov/pubmed/22611049, https://www.ncbi.nlm.nih.gov/pubmed/22735619, https://www.ncbi.nlm.nih.gov/pubmed/22815867, https://www.ncbi.nlm.nih.gov/pubmed/23389292, https://www.ncbi.nlm.nih.gov/pubmed/23446634, https://www.ncbi.nlm.nih.gov/pubmed/23792061, https://www.ncbi.nlm.nih.gov/pubmed/23935582, https://www.ncbi.nlm.nih.gov/pubmed/23996192, https://www.ncbi.nlm.nih.gov/pubmed/24033266, https://www.ncbi.nlm.nih.gov/pubmed/24121126, https://www.ncbi.nlm.nih.gov/pubmed/24398642, https://www.ncbi.nlm.nih.gov/pubmed/24755831, https://www.ncbi.nlm.nih.gov/pubmed/24798971, https://www.ncbi.nlm.nih.gov/pubmed/24825461, https://www.ncbi.nlm.nih.gov/pubmed/25071582, https://www.ncbi.nlm.nih.gov/pubmed/25112650, https://www.ncbi.nlm.nih.gov/pubmed/25287020, https://www.ncbi.nlm.nih.gov/pubmed/25741868, https://www.ncbi.nlm.nih.gov/pubmed/26024779, https://www.ncbi.nlm.nih.gov/pubmed/26025379, https://www.ncbi.nlm.nih.gov/pubmed/26159428, https://www.ncbi.nlm.nih.gov/pubmed/26365338, https://www.ncbi.nlm.nih.gov/pubmed/26416403, https://www.ncbi.nlm.nih.gov/pubmed/26597663, https://www.ncbi.nlm.nih.gov/pubmed/26695521, https://www.ncbi.nlm.nih.gov/pubmed/26716680, https://www.ncbi.nlm.nih.gov/pubmed/26833210, https://www.ncbi.nlm.nih.gov/pubmed/27221532, https://www.ncbi.nlm.nih.gov/pubmed/27255824, https://www.ncbi.nlm.nih.gov/pubmed/27280446, https://www.ncbi.nlm.nih.gov/pubmed/27317329, https://www.ncbi.nlm.nih.gov/pubmed/27437086, https://www.ncbi.nlm.nih.gov/pubmed/27547017, https://www.ncbi.nlm.nih.gov/pubmed/27584680, https://www.ncbi.nlm.nih.gov/pubmed/27657935, https://www.ncbi.nlm.nih.gov/pubmed/27661980, https://www.ncbi.nlm.nih.gov/pubmed/28151915, https://www.ncbi.nlm.nih.gov/pubmed/28684612, https://www.ncbi.nlm.nih.gov/pubmed/28820869, https://www.ncbi.nlm.nih.gov/pubmed/29167213, https://www.ncbi.nlm.nih.gov/pubmed/29201641, https://www.ncbi.nlm.nih.gov/pubmed/30409984, https://www.ncbi.nlm.nih.gov/pubmed/30536387, https://www.ncbi.nlm.nih.gov/pubmed/30798813, https://www.ncbi.nlm.nih.gov/pubmed/31220083, https://www.ncbi.nlm.nih.gov/pubmed/31226389, https://www.ncbi.nlm.nih.gov/pubmed/31469255, https://www.ncbi.nlm.nih.gov/pubmed/31970519
             Polymorphisms: c.118=, c.187=, c.118C>G, c.187C>G, g.8671=, g.8671C>G, g.26090951=, g.26091179=, g.26090951C>G, g.26091179C>G, p.His40=, p.His63=, p.His40Asp, p.His63Asp, n.84=, n.281=, n.309=, n.84G>C, n.281C>G, n.309C>G
     Clinical Significance: RCV000000026.18, RCV000000027.5, RCV000175607.6, RCV000394716.6, RCV000763144.1, RCV000844708.1, RCV000991133.1
                      HGVS:

This question was also asked on Reddit

"What I want to do is to map bases to SNPs (or a nucleotide identifier like NC_000006.11)." - can you please ask this as a separate question; it's inconsistent with the title of this question.

These lines (with no REF or POS) are unmapped reads. If those lines are the entirety of the BAM file, then it is simply another file format to store the raw sequenced reads (as an alternative to FASTQ files). It's somewhat common for BAM files to include all reads, with mapped reads sorted by chromosome and position at the start, and unmapped reads appearing at the end of the file. This means that only the BAM file needs to be transferred and stored, rather than a BAM file and a separate 'unmapped' file (which can change size/content depending on the reference genome and mapper used).
FWIW, the lines of a BAM file can be viewed using samtools view; there's no need to convert to SAM first:
samtools view -h input.bam | less

If you only want to look at the unmapped reads, then you can filter for inclusion using the read flag 4:
samtools view -f 0x04 -h input.bam | less

If you want to exclude unmapped reads, you can filter that flag for exclusion (capital F instead of lower-case f):
samtools view -F 0x04 -h input.bam | less

Its a simple and very clear explanation.

This is an unmapped .bam file.  Someone just converted the fastq to bam format, but no one has aligned these to anything.  So you need to find an aligner program and a genome and align them.


text:
Perform protein structure-based sequence alignment in Python
I am looking for a Python package that performs pairwise structural alignment of protein structures (i.e., PDB files) and returns a sequence alignment. PyMOL is able to do this through the GUI, for example:

For two protein PDBs, one can be aligned to the other with the "align" menu or align command:

With it, the sequences appear aligned in the sequence viewer according to the structural alignment, and an "alignment" object is created, which can be saved with the command:
save A_aligned_with_B.aln, aln_structure_1_to_structure_0
which produces the A_aligned_with_B.aln alignment file in Clustal format (nevermind the ?s, which correspond to special HIS resnames):
CLUSTAL

structure_0  YIYQDSISLPWKVLLVMLLALITLATTLSNAFVIATVYRTRKL?TPANYLIASLAVTDLLVSI
structure_1  EEQGNKL?WAALLILMVIIPTIG-----GNTLVILAVSLEKKLQYATNYFLMSLAVADLLVGL
                                         .*..**..*....**....**...****.****..

structure_0  LVMPISTMYTVT-GRWTLGQVVCDFWLSSDITCCTASIL?LCVIALDRYWAITDAVEYSAKRT
structure_1  FVMPIALLTIMFEAMWPLPLVLCPAWLFLDVLFSTASIM?LCAISVDRYIAIKKPIQANQYNS
             .****......  ..*.*..*.*..**..*....****.***.*..***.**........ ..

structure_0  PKRAAVMIALVWVFSISISLPP-FFWRQAKAEEEVSE-CVVNTD?IL-YTVYSTVGAFYFPTL
structure_1  RATAFIKITVVWLISIGIAIPVPIKGIETDVDNPNNITCVLTKERFGDFMLFGSLAAFFTPLA
             ...*...*..**..**.*..*  .              **....    ........**..*..

structure_0  LLIALYGRIYVEARSRILKQTPNRTGKLLEKKKLMAAR--ERKATKTLGIILGAFIVCWLPFF
structure_1  IMIVTYFLTI?ALQKKAYLVKNKPPQRLTIGKKSVQTISNEQRASKVLGIVFFLFLLMWCPFF
             ..*..*....                              * .*.*.***....*...*.***

structure_0  IISLVMPICKDACW--F?LAIFDFFTWLGYLNSLINPIIYTMSNEDFKQAF?KLIRFK--
structure_1  ITNITLVLCDS-CNQTTLQMLLEIFVWIGYVSSGVNPLVYTLFNKTFRDAFGRYITCNYR
             *.......*   *   .  .....*.*.**..*..**..**..*..*..**

I need a Python package that is easily installable/distributable that accepts two PDBs or protein structures and returns the structural alignment of the sequences somehow. All within a Python script, and not creating intermediate files, e.g. for the alignment.
There are plenty of Python tools that allow structural superimposition and RMSD calculation (e.g., Biopython), but I haven't found any that are able to return a sequence alignment from the structural alignment.
Edit:
Possible answer, thanks to the answer and comments from Matteo Ferla.
import pymol2
with pymol2.PyMOL() as pymol:
    pymol.cmd.load('foo.pdb', 'foo')
    pymol.cmd.load('bar.pdb', 'bar')
    pymol.cmd.align('foo', 'bar', object='aln')
    aln = pymol2.cmd2.pymol.exporting.get_alnstr('aln', _self=pymol.cmd) # I don't know if this line could be cleaner
    # aln is itself a string

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.

Please provide more details to help us in answering your question. For example, showing the alignment that you would like to see from Chimera or PyMOL. Showing the first few lines of any input files is also helpful (or, at least, providing two demonstrative PDB IDs that you are interested in).

This was asked before here —I know because I answered it. You want the primary sequence alignment based on PyMOL's superimposition of tertiary structure.

Ignore my comment. I saw your Q was just updated...

"align" performs a sequence alignment followed by a structural
superposition, and then carries out zero or more cycles of
refinement in order to reject structural outliers found during
the fit. "align" does a good job on proteins with decent sequence
similarity (identity >30%). For comparing proteins with lower
sequence identity, the "super" and "cealign" commands perform
better.

cealign aligns two proteins using the CE algorithm. It is very robust for proteins with little to no sequence similarity (twilight zone). For proteins with decent structural similarity, the super command is preferred and with decent sequence similarity, the align command is preferred, because these commands are much faster than cealign.

yoiu can save a clustal object as a result of cealign, cealign target, mobile [, target_state [, mobile_state
[, quiet [, guide [, d0 [, d1 [, window [, gap_max
[, transform [, object ]]]]]]]]]]   --> object = string: name of alignment object to create {default: (no alignment object)}  https://pymolwiki.org/index.php/Cealign

PyMOL in itself can be used as python package.
The GUI app has its own Python, but in a conda environment you can install it as a package.
Any PyMOL app command can be used as a cmd command. In PyMOL Wiki you see both. They follow a similar pattern and actually have more functionality.
import pymol2
with pymol2.PYMOL() as pymol:
    pymol.cmd.load('foo.pdb', 'foo')
    pymol.cmd.load('bar.pdb', 'bar')
    pymol.cmd.align('foo', 'bar')
    pymol.cmd.save('foobar.aln')

For more discussion of how the alignment output object is used see Assessing PyMol sequence alignment object

Hi! Thanks so much for your answer. Sorry it took me so long to say anything, I was testing all the options. Also, I'm sorry I don't have enough reputation to upvote your answer. Is there any way I could save the alignment in a Python text stream, as in an io.StringIO() stream? So that I don't save any file in the process. I made some attempts but I couldn't hack it and I didn't see anything related in your previous answer that you linked. Thanks in advance!

Nah, the important thing is that marked as solved, thus helping the healthiness of the site, which struggled to leave beta because of high unsolved cases —so thanks!
I don't know off the top of my head, but https://docs.python.org/3/library/tempfile.html may be of use. There is likely a function that is not imported by the top of cmd module, but is present in one of the submodules. Maybe one of the searches in the open-source PyMOL repo might help https://github.com/schrodinger/pymol-open-source/search?q=.aln (sorry I'm tight for time atm)

Thanks for the pointers! I finally found out how to retrieve the alignment string, I'm going to edit my question to put the final answer in it.

@MatteoFerla any way to use it in python without anaconda enviroment ? https://stackoverflow.com/questions/63958085/how-to-embed-the-pymol-module-into-python

Hi @pippo1980 maybe a separate question? Generally not advised, you could use venv and pip, but historically there were compatibility issues between conda and pip (probably resolved)

Anyway thanks @FranchoNerínFonz for finalising the post. Very interesting (and useful) question and answer indeed.

@pippo1980 The classic solution was installing the OS pymol —which is not installed via conda but with wizardry and  the deprecated python setup.py install (not pip install .) https://github.com/schrodinger/pymol-open-source/blob/master/INSTALL —alternatively, there's always the nightmarish über-hack of copypasting of site-packages folders from a conda installation to a vanilla Python... But that's a big no no due to irreproducibility.


text:
Calculating the effect of differences between two groups of nucleotide sequences at every position
I'm investigating two different groups of short nucleotide sequences. The behaviour of a certain DNA-binding protein is likely influenced by common differences between these two groups leading high or low activity. Some positions are hypothesized to be more important than others. All sequences have the same length and are very similar.
I'm trying to find out how well changes at certain sequence positions explain the between-group difference and which positions are of the highest importance for this relationship. By looking at the sequence logo for both groups I can easily identify crucial positions. However, I'm at a loss how to actually score the effect each position has and how well it correlates to membership in either group. I thought about calculating Cramers V at each position but that doesn't seem to be very useful.
Can someone here point me in the right direction?

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.

are you talking about a sequence alignment of homologous sequences? Can you include a sample of what your data look like? Recall that the logo is based on sequence entropy, why not use that measurement?


text:
Retrieve protein sequence from Mgnify given only accession code
I only have the accession codes of several proteins from the MGnify database (https://www.ebi.ac.uk/metagenomics/). I would like to retrieve the full amino acid sequence data from the database, but I have not been able to find a way to do this using the web server or available API.
The website provides sequence search with the HMMER tool (https://www.ebi.ac.uk/metagenomics/sequence-search/search/phmmer). However, I need the sequence data to search with HMMER.
How can I get the protein sequence only with the accession code?
P.S: An example of the accession code is MGYP001242958860.


text:
Clustering individuals by gene presence/absence
I have a binary matrix with individuals as row names, and gene name as column name.
So, if the gene is present in an individual, we have 1, otherwise 0.
I would like to cluster individuals based on the gene presence/absence pattern.
Is it right to perform an MCA? Or a logistic PCA is needed?
Thanks,
Marco

Would you consider hierarchical clustering, i.e. making a tree? However you need to present a test data set and I could give you the tree with an explanation. PCA is okay, but trees have advantages, for reasons I could explain.

I believe the answers to your question are here:
https://help.ezbiocloud.net/gene-presenceabsence-analysis/
https://www.biostars.org/p/66355/


text:
Heatmap of motifs present in protein sequences
I have a file containing multiple protein amino acid sequences. I am interested in visually presenting the location of amino acid motifs PQG, QQG and AQG in the sequence in the form of a heatmap in R. Every location of the motifs should be represented as a differently coloured barcode. Example sequence follows,
>cat archae.fasta
>A0A075FKY9
MRRITLPKGIEISGKSPFGTSSGVYEVVISISENKPSDASIASNSFNSLWKDNFHLRVSG
DHFSQVLGSGSNPIPDNVFDLGSVWIVIQDQFSPVHTSFQFNISPSNQPASKPEEKPRRV
ETSTVKRTITKPGRLGPPGGKGPRGPSGYPGIKGDKGGRGPTGDKGDKGDKGIPGPQGEK
GKTGQAGDKGDKGITGPPGEKGPRGSTGPPGDKGDKGVQGKIGDKGLTGTTGPVGDKGTQ
GPQGPPGGKGLTGIPGPQGDKGTKGPLGPVGERGPTGKPGESGLQGPQGIQGAQGERGHS
GPPGPQGEKGLIGDKGEIGERGTRGPPGPPGEKGSQGGMSDEGKRLIKELLELLASKNII
STEEQIKLTSYLY
>A0A075FMB0
MPGDKGTSGIRGVPGDKGDKGPQGPPGDKGLTGPVGVPGEKGPTGTSGVQGEKGIQGTPG
PLGEKGVTGPAGDKGPIGPPGPAGSKGITGPAGPLGDKGTIGPPGPLGDKGSKGPEGPTG
DKGPQGPQGPAGAKGLTGVPGPQGEKGEKGPLGPIGEKGTTGATGPPGDKGPQGPQGTQG
ERGTTGPSGEKGPQGPQGIQGPQGERGPTGPIGSIGETGPIGGQGPVGPAGPRGPPGPPG
EKGPSGGMSSEQKAIFKELLEILTTKEIISTEEQIKLMSYLY
>A0A075G2F8
MPKGIEISGKSPFGTSSGVYEVKISVSQSKPSNDSITNNSFNSLWKDNFHLRVSEGYFSE
ILGSDSNPIPDNVFDLGSVWIVIQDQFSPVHTSFEFKISQDSKSIPTPKEKPRRVDTSSV
KKVRTIPTRQGPTGSKGIPGERGYPGTKGEKGGKGPTGDKGDKGDKGVPGPQGEKGKTGS
TGDKGDKGITGSPGEKGPRGPTGPPGDKGDKGVQGKIGDKGLTGTTGPIGDKGPQGPQGP
PGDKGLTGIPGPQGDKGGRGTIGPSGERGPTGRPGESGLQGPQGIQGPQGERGHSGPPGP
QGEKGIIGDKGEIGERGPRGPPGPPGEKGSQGGMSDDGKILFKELLQLLASKNIISTEEQ
IKLTSYLY
>A0A075GLR7
MPKGIEISGKSPFGTSSGVYEVRISISQDKPSDDGITNNSFNSLWKDNFHLRVSEGHFSE
ILGTDSNPIPADVFGLGSIWIIIQDQFSPVHTSFELHISQDNKLLSVPQGKPKPVNATKV
QKVRAQPGRPGPQGEKGGKGPAGYPGTRGEKGGRGITGDKGDKGDKGIPGPQGEKGIKGP
AGDKGDKGITGPVGPSGPKGSTGPPGDKGDKGVQGKSGDKGLTGVTGPLGDKGTQGPQGP
PGPKGLTGVPGPQGDKGMRGSIGPAGERGPTGKPGDPGLPGPQGIQGPQGERGHTGPQGP
PGEKGLVGDKGEIGDRGPRGPPGPPGEKGSQGGMSEESKRLVKELLELLASKNIISTEEQ
IKLTSYLY

As an example, please see a two-colour heatmap. I need something similar but with three colours corresponding to each of the three motifs.

I looked at the heatmaply package but the biggest problem is that my matrix is not numeric. Is there way to make heatmap using text strings in R?
Motif is defined as any three character string that matches the general pattern either PQG, QQG and AQG in the uploaded sequence. y-axis is the fasta code appended after > (A0A075FKY9 for example). x-axis would depend on how the data is converted into a matrix.

Firstly, the problem with your data is that it ain't aligned. Without that your motif positions are not relative from one protein to another. Basically you need to align it and there's no wriggle room.
Secondly, its easy you simply parse the motifs you want and their location within the string - after it is aligned. In Python you would use finditer within a comprehension. This will go into a dataframe and subject to a heatmap package, I'm not a 'heatmap' person but it would be easy to graph where the x axis is the position where each motif is found and a qualitative score given for each type of motif (1-3 here). I don't know the equivalent in R.
There is an issue whether an indel will remove a motif and you would need to check for that by running the same command against an unaligned protein and seeing if the frequency matches.
My advice is to ask a new question with some code. I dunno whether someone has coded this already, - its easy to do. Personally, I don't think there's sufficient requirement, but I could be wrong.


text:
How to remove the unpaired reads in sam/bam files?
I have sam and bam files for the chimeric reads, which come from two different parts of the genome (For example, the first half of the read from part of Chromosome 1 and the second half of the read from part of Chromosome 3). I have removed the low mapping quality reads (MAPQ <= 30) and reads with high mismatches. This might result in some unpaired reads which need to be removed before further analysis. May I ask whether there is any way to remove the unpaired reads in sam/bam files? Thanks!

The -f option of samtools view is for flags and can be used to filter reads in bam/sam file matching certain criteria such as properly paired reads (0x2) :
samtools view -f 0x2 -b in.bam > out.bam

Thanks. If I remove the low mapping quality reads and reads with high mismatches after the alignment, will the flags be updated? If they are not updated, there will possibly be unpaired reads even if I use flags for filtering in my point of view.

That is a good question and there is a good change those flags won't be updated, it would need some testing. Here is a post on seqanswers that might help you.


text:
How does split reads look like in sam files?
I used bwa mem to align the DNA with the reference genome. If there are split reads (chimeric reads, come from two different parts of the genome), will it be split into two lines rather than one line?

Split reads as in spliced or split reads as in a read with a supplemental alignment? "Split reads" has many possible meanings.

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.

split reads here refer to chimeric reads which come from two different parts of the genome. For example, the first 1-70 bp is from a region in chromosome 1, 71-150bp is from a region in chromosome 3.

I'm assuming you mean supplementary/chimeric alignment? The SAM Format Specification has a really detailed explanation as well as its Optional Fields Specification.
In general though, you will have an additional line in the SAM file to represent a chimeric alignment yes. For a single read, one line will represent the primary alignment; this is arbitrarily chosen. The other line will represent the split/supplementary alignment and have its supplementary alignment flag set. If there are multiple chimeric alignments, each one will have its own line.
You check if a primary alignment has any supplementary alignments associated to it by checking the SA tag which will list the location(s) of any supplementary alignment(s).


text:
Can I Incorporate svaseq() into GSEA/GSVA analysis?
I understand GSEA/GSVA can take microarray-like expression matrix such as the output of voom() or vst(). However, I have a question on how we can also use svaseq() variables, correct the output from voom() or vst(), and apply GSVA on it.
I assume this is possible since I have read papers showing PCA after correction by SVA (such as this one and this one in their supplementary figure), but my statistical and coding knowledge does not allow me to recreate what they do. Is there a known pipeline that can do this
Edit: to clarify, I would like to generate a normalised matrix (e.g. rlog, or out put of voom from limma, or vst from DESeq2) that is adjusted for batch effect. ComBat, for example, can provide what I need. However, since I am doing my DESeq2 DGE analysis with svaseq(), I want to adjust it with svaseq variables. I am wondering if there are pipelines that can use svaseq vairables (or any varaibles) and adjust the value in the normalised expression matrix.

Use removeBatchEffects from limma. The input counts should be on log scale, so vst and rlog are ok. The covariate argument can take the surrogate matrix from svaseq.
https://rdrr.io/bioc/limma/man/removeBatchEffect.html


text:
How should I interpret DGE results if only one HLA-A gene shows up as significant but not the others?
I have done a DGE recently and have been looking at the DGE list. One of the genes is HLA-A. However, when I dug deeper I realised there are hundreds of HLA-A genes with unique ENSEMBL number (of course it's HLA...). A few other HLA-A genes also show up on the table but are not significant at all.
Considering the way Salmon deals with multimapping reads, can I assume these expressing genes are different enough that they shouldn't be consider the same? Or should the expression of all HLA-A genes considered as one single gene for the purpose of identifying differentially expressed genes?
I do understand that this may be more of a biological question than bioinformatics, but I also want to see if it is a practice pouring the HLA genes (or other variable genes/loci of similar nature) with different sequences together when doing DGE, and whether the statistics would be valid if I do so?

It is bioinformatics IMO because HLA is a massive area of human population genetics and its complex. As a personal opinion I think you need further supporting evidence.
If I understand the scenario correctly,

you are confident the data is correct (see below),
a given HLA-A genotype is significant above other HLA-A and HLA-B genotypes

Thats publishable result, cutting-edge potential, you just need to find the given genotype and represent the association within a given disease association study. Yes your result is absolutely of interest particularly if you consider the wider output of association studies (HLA-B). It is trendy to identify a given HLA to a given association. If this is non-pathogen related they HLA exhibits very strong population skews.
I have put a large caveat below, because:

you shouldn't underestimate the complexity of this research field and
I don't really understand the specifics of your study.

If you want to ditch or include HLA it would require carefully defined biological aims. It is a specialist area of investigation IMO (below).
HLA diversity and its link to population shifts in infection susceptibility is subject to repeated studies. Different populations show very distinct allelic frequencies. For example, "HLA-B*40 Allele Plays a Role in the Development of Acute Leukemia". Distinct link between African population and malaria .. These just being a very small sample in a hugely active area of investigation. The idea is that specific genotypes have very specific outcomes in human health.
Rationale The single gist is HLA this represents the cellular immune response and is considered essential for infection control, with the notable exception of pathogenic bacteria. Each genotype is considered to be uniquely functional and again is well known to exhibit population skews. It is beyond 'trendy' is huge and established field of investigation.
Caveat If you were looking for 'wiggle room' to avoid HLA-A diversity and allelic frequency and all the associated work it requires ... or you specifically want to make it a standout result ... then the angles that might work are ...

alot of examples of association focus more on HLA-B than HLA-A. This is a rule of thumb and might be disputed. So if you're going to ignore HLA-A you could perform an extensive literature review to see the extent that associations of HLA-B predominate. I think you might find that holds. So its either really cool cause its an outlier, or you can ditch it because its an outlier.
If pathogenic bacteria were the focus of your investigation, then most would accept HLA is not relevant.

I'm sure there's other angles of justification, e.g. "I recognise HLA-A .... but this will be subject to further investigation by (putative) collaboration".
Advice If I was tackling this data I would contact a bioinformatics immunologist, because they will have databases and pipelines dedicated to the exact task you require. Personally, I would be reluctant to take this on without specialist assistance, at a minimum someone with current experience. Essentially, what I'm trying to say, is you will need to understand the HLA field because crossing into it isn't trivial. The situations where I could be wrong are that its simply so complex that no-one really understands it, in which case your good to go.

That's really well explained! Thanks @M__ ! So the research question I have is basically the transcriptomic difference between haematological malignancy which enter remission and those which relapse. This may sound stupid but are these unique ENSEMBL IDs of HLA-A just a different allele/genotype/serotype (I am not familiar with the terminology)? And since each person has only a pair I can tell which alleles they have by RNA-Seq (theoretically)?

I don't know this data set. The data should be individual HLA-A genotypes if correctly annotated, which they may not have done. The second part of the question is definitely yes you can tell if correctly annotated (I don't understand 'one pair'). This may help https://raw.githubusercontent.com/ANHIG/IMGTHLA/Latest/alignments/A_nuc.txt  http://hla.alleles.org/data/hla-a.html  Finally, there are no serotypes in HLA, for an explanation of terminology https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3715123/ .

Thanks a lot! Just to add in case people are reading, the unique ENSEMBL IDs are from MHC contigs from different cell lines. I will try to see if it is possible to get the genotype of the samples by WGS and align RNA-Seq reads accordingly. To clarify, when I say pair I assume like other genes in the genome there are two alleles hence a pair of the genes. I will look up more about HLAs and MHCs.

@Kento MHC  is basically the same thing, so that data looks fine. However, if this is a mouse model then (obviously) they don't have HLA just MHC.

You need specialized tools to study HLA genes. Generic tools won't work well. Before other analysis, you should get the correct genotype first. Try arcasHLA. It has been a popular HLA genotyper for RNA-seq data. When you know the genotype, I guess you may get rid of other HLA-A alleles and redo the differential analysis if you are interested in HLA-A. I haven't done this, though. Just a thought.