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.


text:
change/edit FASTA headers
I have a FASTA file with 97,000 entires arranged like this :
>tr|A0A075B6G3|A0A075B6G3_HUMAN Dystrophin OS=Homo sapiens OX=9606 GN=DMD PE=1 SV=1

I want to have an output like this:
>A0A075B6G3|1|training

I am aware of many one liners in sed, perl and awk but none are working . I have something like this :
awk '/>/{$0=">ID|1|training"}1'

Kindly help.

Where is the |1| coming from? Is that an incrementing value that will change for each line? Is it the PE=1? Or maybe the SV=1?

One way using awk:
awk 'BEGIN { FS=OFS="|" } /^>/ { print ">" $2, "1", "training"; next }1' file.fa

This sets the field separator (FS) and the output field separator (OFS)  to a pipe character. Then for each header line, we can just print the fields we want. The output fields will be delimited with pipe characters since we already set OFS to a pipe character in the BEGIN block. We can use next here to force awk to immediately stop processing the current record and move on to the next record, which is just the next line, since, by default, awk's record separator is a newline character.
Since the 1 on the end just returns true, it performs the default action, which is to print the current record stored in $0. It is just the same as:
awk 'BEGIN { FS=OFS="|" } /^>/ { print ">" $2, "1", "training"; next } { print $0 }'

Which is the same as:
awk 'BEGIN { FS=OFS="|" } /^>/ { print ">" $2, "1", "training"; next } { print }'

Thank you @ Steve.

Nice! You could shorten it to something like awk -F'|' '/>/{printf ">%s|1|training\n",$2; next}1' file.fa or even awk -F'|' '{ (/>/) ? $0=">"$2"|1|training" : $0=$0}1' file.fa if you want.

perl -p -i -e 's/>[a-z]+?\|([A-Z]+[^\|]+?\|).*SV=([0-9]+)/>$1$2\|training/g' myfilename.fa

Not tested yet. Just check the sequence ID starts with a capital A-Z character. Probably wise to test it without the -i

Many thanks @M_

Thanks @user1738234 , I'll test it later. I assume that the 1 refers to the SV number. Is that correct?

Here's another simple perl solution:
$ perl -lpe 's/\|.*(\d+)/|$1|training/ if s/>.+?\|/>/' file.fa
>A0A075B6G3|1|training

Note that you didn't tell us where the 1 came from, so I am taking the last numbers on the line, the SV=1 in your case. If you want to hard-code the 1, you can do so easily with:
perl -lpe 's/\|.*/|1|training/ if s/>.+?\|/>/' file.fa

The -lpe options mean "print each line, adding a newline character to each print call, after applying the script given by -e". So this will read the file line by line, do its thing and print the result.
The s/old/new/ is the substitution operator and will replace old with new. In this case, we are replacing a | character (\|) and everything after it (.*) until the last set of numbers on the line which we capture and save as $1 with (\d+). This substitution will only happen if we have already successfully substituted the leading > and everything until the first | (>.+?\|) with just a >.
If you want to edit the original file in place, use the -i flag (but give it a value, .bak in the example below, since that ensures a backup file with the original data and that extension will be created):
perl -i .bak -lpe 's/\|.*(\d+)/|$1|training/ if s/>.+?\|/>/' file.fa

source: Stack Exchange

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

sort: votes


text:
How to tune and use the MetaVolcanoR package
I conducted a differential expression analysis over several datasets, using LIMMA, each one on its own. For each dataset, I have a data frame of all the genes, with the LFC and the p-value. Those data frames are called deg1, deg2, deg3 ... and so on.
I want to use the MetaVolcanoR package, using the votecount_mv approach since it is the only approach that I can pick the p-value and LFC thresholds that I desire.
The problem is I'm not sure how to tune the parameters. There are several problems that I'm facing and I would really use some help please:
How do I set the log fold change threshold if I want only LFC > 0.5 and LFC < -0.5?
Is it log fold change = 0.5 ? just making sure.
What is the metathr parameter? I read about it and still don't understand what it does. I set it as 0.01 for now.
And the main problem, the plot is showing much more gene number than there actually is. You can see in Dataset6 for example, the plot says it has over 50,000 genes while in reality, it has something like 18,000 genes. The same for Dataset8, why is this happening?
This is the plot:

This is my code:
totalDEG = list(Dataset1 = deg1, Dataset2 = deg2,
                Dataset3 = deg3, Dataset4 = deg4, Dataset5 = deg5, Dataset6 = deg6,
                Dataset7 = deg7, Dataset8 = deg8 ,Dataset9 = deg9, Dataset10 = deg10,
                Dataset11 = deg11, Dataset12 = deg12, Dataset13 = deg13 , Dataset15 = deg15,
                Dataset16 = deg16, Dataset20 = deg20, Dataset21 = deg21)

totalDEG = map(totalDEG, ~ .x %>% rownames_to_column("symbol") %>% `row.names<-`(.$symbol))

meta_degs_vote <- votecount_mv(diffexp=totalDEG,
                               pcriteria='P.Value',
                               foldchangecol='logFC',
                               genenamecol='symbol',
                               geneidcol=NULL,
                               pvalue = 0.05,
                               foldchange = 0.5,
                               metathr=0.01,
                               collaps=FALSE,
                               jobname="MetaVolcano",
                               outputfolder=".",
                               draw='HTML')

head(meta_degs_vote@metaresult, 50)
meta_degs_vote@degfreq

This is a small sample of the deg13 dataframe, just as an example of how it looks:
structure(list(logFC = c(3.39094029649578, 2.91185452846888,
-4.18752593306854, 4.15866027194434, -3.46058235630518, -2.87208933878881,
-3.08172118179278, -2.7458737595494, -4.086372553458, -2.68889593085408,
-3.03835320252487, 2.61330670793883, -2.48123655627091, 1.49263378072949,
1.95754205818365, -2.70230318294949, -2.9154447157702, -3.38110311321268,
-2.35284609598761, -3.61327208834916, 2.03692346164301, -1.70931971282814,
2.27925549238084, -3.79524834465106, 1.64515339181668, -4.30959811283984,
3.57930610863801, -2.54394115096987, -1.98143797060853, 4.03178891878422
), AveExpr = c(6.4934411845873, 6.79507589173635, 3.87626322696796,
5.02990511836996, 3.37154863566726, 4.54352970704667, 2.51617183960423,
5.36760482831628, 8.74162861573995, 3.52708473144296, 1.56671205994866,
5.87684931921495, 4.44218678637856, 5.12219572880673, 7.17539095745285,
3.53733757517713, 3.35148844504943, 5.15160050823647, 4.75910049955374,
2.76067393839851, 7.59404366884022, 3.26758602409262, 6.6742932845645,
2.02865432630183, 7.35343444625846, 1.63566796485655, 6.63046957559406,
4.0296762052982, 5.52320581646106, 5.21519413005468), t = c(7.67283765970783,
7.48874311179818, -6.82705440075341, 6.17183091062419, -6.04487485029223,
-5.8487763400982, -5.92000790680834, -5.64958317461715, -5.61111489307273,
-5.58440722941211, -5.78651741270403, 5.34175048353288, -5.31088208376663,
5.08100365503679, 5.03973023175764, -5.07773225272668, -5.04632997176375,
-4.96736474738699, -4.82583305184561, -4.92525789154973, 4.78892270501915,
-4.75294466621978, 4.64502680172295, -4.88600418339975, 4.62243262623093,
-4.78866983562446, 4.56376665911121, -4.57063797766699, -4.48421852802591,
4.47493043300733), P.Value = c(8.61557824744483e-06, 1.08572492426076e-05,
2.57808593137873e-05, 6.41321288302935e-05, 7.70283126846018e-05,
0.000102673342223582, 9.24388613776644e-05, 0.000138234657017388,
0.000146500473045098, 0.000152546919410067, 0.000112607650587545,
0.000221306832937214, 0.000232172444273951, 0.000333174926637586,
0.000355778600664596, 0.000334909952712188, 0.000352058129044564,
0.000399405716020234, 0.000501862041539822, 0.000427357183393412,
0.000532900955175219, 0.000565105101377844, 0.000674574457630558,
0.000455275577064366, 0.000700192843853878, 0.000533120436605821,
0.000771595292301189, 0.000762850982694916, 0.00088083209108589,
0.000894604612049976), adj.P.Val = c(0.105689892752164, 0.105689892752164,
0.167309183326708, 0.24749466449955, 0.24749466449955, 0.24749466449955,
0.24749466449955, 0.24749466449955, 0.24749466449955, 0.24749466449955,
0.24749466449955, 0.322868951254968, 0.322868951254968, 0.364560714544159,
0.364560714544159, 0.364560714544159, 0.364560714544159, 0.388801494259897,
0.397198753468022, 0.396200809689826, 0.397198753468022, 0.397198753468022,
0.397198753468022, 0.397198753468022, 0.397198753468022, 0.397198753468022,
0.397198753468022, 0.397198753468022, 0.397198753468022, 0.397198753468022
), B = c(3.51320342275366, 3.35460378796296, 2.44903607778273,
1.62363531926393, 1.57839963166078, 1.44388223679157, 1.36898233186455,
1.24528027435772, 1.14177006290097, 1.04403303682497, 0.945243757754988,
0.859404654981856, 0.797031901979986, 0.496397382158138, 0.450148400219772,
0.418557695883232, 0.3417613688148, 0.324284286390148, 0.136478287041671,
0.0982449107580958, 0.0925437909507512, -0.0280517206307875,
-0.112494849608916, -0.129424832005619, -0.147717311976667, -0.200056933894112,
-0.234222113297706, -0.244412573722838, -0.351160223852428, -0.386270016806261
)), row.names = c("FOS", "DUSP1", "TP53AIP1", "FOSB", "TFAP2C",
"KIF26A", "GRHL3", "CTAGE4", "FCGBP", "FBXO27", "AC087289.3",
"SH3BP5", "GABRE", "ADAMTS4", "DDHD1", "KCNS3", "ARHGEF16", "SERINC2",
"OPLAH", "WNT7B", "GSN", "OVOL2", "CD163", "SVOPL", "ZFP36",
"SPESP1", "FGD5", "FAAH", "ERBB3", "PARM1"), class = "data.frame")

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


text:
How to select the best ligands in a virtual screening matrix
I have the results of a virtual screening experiment using docking simulations with Autodock Vina. The result is a matrix of 7 (proteins) by 28000 (ligands) with the calculated binding energies for each protein/ligand pair:

Ligands - Proteins
Protein1
Protein2
Protein3
Protein4
Protein5
Protein6
Protein7

Ligand1
-5.4
-5.9
-5.0
-5.2
-5.8
-5.7
-4.5

Ligand2
-9.5
-8.6
-8.9
-9.0
-9.2
-7.1
-6.4

Ligand3
-9.5
-8.3
-8.2
-9.9
nan
-6.8
nan

Ligand4
-9.4
-7.8
-8.0
-6.9
-9.5
-1.3
-7.8

Ligand5
-9.3
-6.9
-8.5
-7.9
-9.2
-6.8
-6.9

---

Ligand28000
-6.3
-7.9
-5.4
-2.7
-10.1
-9.4
-3.9

nan results are either failed dockings, or positive values.
The scores are not only affected by how well a ligand docks with each protein, but also bigger (higher molecular weight) ligands tend to show better (more negative) scores. However, this is in part an artifact and it does not imply that bigger molecules are better and smaller ones are worse.

Typically, for a virtual screening with just one protein, you would select the top ranking ligands and that would be it (e.g. keep the top 1% and test them). However, this results in some ligands selected for most proteins because they invariably score better due to their size.
In this case, I could also take into account the variability of each ligand because I have observations for different unrelated proteins (i.e. see if a particular ligand systematically binds with better scores across proteins).
Is there any way to take into consideration both factors when selecting the best ligands?
I thought of calculating Z-scores 1) across rows and 2) across columns, and select the ligands for each protein that have less than -1.65 standard deviations (p=0.05) and that also fill the same criteria across proteins. However, there are only 7 observations for each ligand to calculate the Z-scores, and it does not feel correct to give it the same weight to the Z-scores calculated from the 28000 observations for each protein. Plus, in the end, protein docking effect is more important than the ligand effect.
Any suggestions would be greatly appreciated!
Ligand efficiency looks like a promising suggestion. However, upon trying (calculating ∆G/MW or ∆G/n.heavy_atoms) resulted in very small ligands with mediocre or bad energies going to the top (e.g. urea with -3.1 became a Z-score of -8.5 SDs).

Urea?! —what dataset are you docking and what's the application? I had assumed it was Enamine REAL or Wuxi biologics dataset for purchasing for drug discovery.

If it is to find the native ligand of an enzyme, the task is completely different, but the monetary aspect is worse —say the native ligand is something like phosphoribosyl pyrophosphate and you have 10 candidates that are close, that will cost as much as the GDP of a small state. Therefore ancillary checks are needed like operon guilty-by-association etc.

I seem to have glitched the system by adding an extra comment and my previous comment has shrunk. If it's for drug discovery I'll re type it... :shrung:

Sorry Matteo I don't know where my comment went either! I used the ZINC database of natural compounds here, with the aim of drug discovery. It is expected to find things like urea I guess, and of course any result like this would be manually curated

My comment was about cost and success rate: smaller compounds are cheaper and more likely to be true so LE + a lower cutoff of 150–200 Da depending on the assay or co-crystal soak. So generally smaller compound hits are elaborated upon until one has a lead and jumping to the final product is a waste of cash. ZINC natural products is not make-on-demand & may require a chemist to source the plant (read: import paperwork) and extract the product (read: low yield). Plus they may not be not "drug-like" (cf. Lipinski's Ro5 vs. nature's proclivity towards PO4, CoA or glycosyl groups as a handle).

What I did was to keep any ZINC "natural-products" entry also described as "for-sale", which (in theory...) should be available for buying from distributors, and also discard any that violated 2 or more Lipinski's rules. But yeah, still it left with around 28000 compounds. I will try to see if I can use LE without boosting the "very" small compounds.

Ah. 28k is more than I thought. I find natural compounds awesome to stare at to see their biosynthesis, but I'd worry the routes for expansions may be more isolated than the modular building block approach —say you have a cyclised sesquiterpene (farnesol (3x isoprene) is 222 Da) and you would like to change a substituent it would need to exist in nature or on sale or easily modifiable (Codexis & co. sell plates of exciting enzymes but AFIK that is for production optimisation not exploration).


text:
how to change the pseudo time trajectory direction?
Monocle is a comprehensive bioinformatics package which has the function to run the pseudo-time trajectory and changes in direction.
I wish to replicate this analysis. Like in the figure I would like the right side to be the pseudostar.
Thank you so much for your guys help

Thanks for the figure. Its a standard PCA with time described by a colour gradient. IMO it looks fairly easy to do outside Monocle, just my opinion though.


text:
How does one convert "flowJO" flow-cytometry cell-count, etc, bead, etc data into .json or relational-database-management-system compatible files?
How does one convert "flowJO" flow-cytometry cell-count data into ".json" or relational-database-management-system compatible files?
The way flow-jo exports is in wide rather than long format. It also exports everything per subject/animal/lab-rat per specific timepoints rather than aggregating. I need the aggregated files to plot groups.
It's also not shaped in a "primary-key", column 1, column 2, column 3, or .json (NO_SQL) format.
Running statistical programs or databases is intractable with the export options Flowjo gives.
I am looking for some way to convert it, or maybe the problem lies in flowjo. It does not support the export formats I need.
I amk open to other flow-cytometry software that offers better integration with no-sql/sql database.

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.


text:
Do proper .fcs flow-cytometry files specify the gating hierarchy, so one may calculate cell-counts without specifying it themselves?
Do proper .fcs flow-cytometry files already have the gating hierarchy, so one may calculate cell-counts without extra work?
https://github.com/whitews/FlowKit/
# setup some file paths for our data
base_dir = "data/8_color_data_set"

sample_path = os.path.join(base_dir, "fcs_files")
wsp_path = os.path.join(base_dir, "8_color_ICS.wsp")
session = fk.Session(sample_path)
session.import_flowjo_workspace(wsp_path)
# that's it, let's see what sample groups were loaded
session.get_sample_groups()
['default', 'All Samples', 'DEN']
# I happen to know all the "real" analysis is within the "DEN" group
sample_group = 'DEN'
# and the sample that are included in the group
session.get_group_sample_ids(sample_group)
['101_DEN084Y5_15_E03_009_clean.fcs',
 '101_DEN084Y5_15_E01_008_clean.fcs',
 '101_DEN084Y5_15_E05_010_clean.fcs']
# so what is the gating hierarchy for this group
print(session.get_gate_hierarchy(sample_group))
root
╰── Time
    ╰── Singlets
        ╰── aAmine-
            ╰── CD3+
                ├── CD4+
                │   ├── CD107a+
                │   ├── IFNg+
                │   ├── IL2+
                │   ╰── TNFa+
                ╰── CD8+
                    ├── CD107a+
                    ├── IFNg+
                    ├── IL2+
                    ╰── TNFa+
# looks good, let's go ahead and analyze them (using verbose mode to see each gate as it's processed)
session.analyze_samples(sample_group, verbose=True)
#### Processing gates for 3 samples (multiprocessing is enabled - 3 cpus) ####
101_DEN084Y5_15_E03_009_clean.fcs: processing gate Time
101_DEN084Y5_15_E03_009_clean.fcs: processing gate Singlets
101_DEN084Y5_15_E03_009_clean.fcs: processing gate aAmine-
101_DEN084Y5_15_E01_008_clean.fcs: processing gate Time
101_DEN084Y5_15_E01_008_clean.fcs: processing gate Singlets
101_DEN084Y5_15_E01_008_clean.fcs: processing gate aAmine-
101_DEN084Y5_15_E05_010_clean.fcs: processing gate Time
101_DEN084Y5_15_E03_009_clean.fcs: processing gate CD3+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate Singlets
101_DEN084Y5_15_E05_010_clean.fcs: processing gate aAmine-
101_DEN084Y5_15_E03_009_clean.fcs: processing gate CD4+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate CD3+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate CD4+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate CD3+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate CD8+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate CD4+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate CD8+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate CD107a+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate CD8+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate CD107a+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate IFNg+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate CD107a+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate IFNg+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate IL2+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate IFNg+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate IL2+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate TNFa+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate IL2+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate TNFa+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate CD107a+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate TNFa+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate CD107a+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate IFNg+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate CD107a+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate IFNg+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate IL2+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate IFNg+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate IL2+
101_DEN084Y5_15_E03_009_clean.fcs: processing gate TNFa+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate IL2+
101_DEN084Y5_15_E01_008_clean.fcs: processing gate TNFa+
101_DEN084Y5_15_E05_010_clean.fcs: processing gate TNFa+
# and a look a the results
session.get_group_report(sample_group)
sample  gate_path   gate_name   gate_type   quadrant_parent parent  count   absolute_percent    relative_percent    level
0   101_DEN084Y5_15_E03_009_clean.fcs   (root,) Time    RectangleGate   None    None    283968  99.999648   99.999648   1
1   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time)    Singlets    PolygonGate None    Time    236780  83.382341   83.382635   2
2   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets)  aAmine- PolygonGate None    Singlets    161823  56.986150   68.343188   3
3   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-) CD3+    PolygonGate None    aAmine- 132200  46.554377   81.694197   4
4   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+)   CD4+    PolygonGate None    CD3+    81855   28.825330   61.917549   5
5   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+)   CD8+    PolygonGate None    CD3+    46965   16.538777   35.525719   5
6   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) CD107a+ RectangleGate   None    CD4+    59  0.020777    0.072079    6
10  101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) CD107a+ RectangleGate   None    CD8+    814 0.286651    1.733206    6
7   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) IFNg+   RectangleGate   None    CD4+    598 0.210586    0.730560    6
11  101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) IFNg+   RectangleGate   None    CD8+    1224    0.431033    2.606196    6
8   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) IL2+    RectangleGate   None    CD4+    64  0.022538    0.078187    6
12  101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) IL2+    RectangleGate   None    CD8+    127 0.044723    0.270414    6
9   101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) TNFa+   RectangleGate   None    CD4+    288 0.101420    0.351842    6
13  101_DEN084Y5_15_E03_009_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) TNFa+   RectangleGate   None    CD8+    821 0.289116    1.748110    6
0   101_DEN084Y5_15_E01_008_clean.fcs   (root,) Time    RectangleGate   None    None    290166  99.997932   99.997932   1
1   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time)    Singlets    PolygonGate None    Time    239001  82.365287   82.366990   2
2   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets)  aAmine- PolygonGate None    Singlets    164655  56.743931   68.893017   3
3   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-) CD3+    PolygonGate None    aAmine- 133670  46.065782   81.181865   4
4   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+)   CD4+    PolygonGate None    CD3+    82484   28.425899   61.707189   5
5   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+)   CD8+    PolygonGate None    CD3+    47165   16.254153   35.284656   5
6   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) CD107a+ RectangleGate   None    CD4+    68  0.023434    0.082440    6
10  101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) CD107a+ RectangleGate   None    CD8+    73  0.025157    0.154776    6
7   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) IFNg+   RectangleGate   None    CD4+    4   0.001378    0.004849    6
11  101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) IFNg+   RectangleGate   None    CD8+    2   0.000689    0.004240    6
8   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) IL2+    RectangleGate   None    CD4+    6   0.002068    0.007274    6
12  101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) IL2+    RectangleGate   None    CD8+    10  0.003446    0.021202    6
9   101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) TNFa+   RectangleGate   None    CD4+    21  0.007237    0.025459    6
13  101_DEN084Y5_15_E01_008_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) TNFa+   RectangleGate   None    CD8+    8   0.002757    0.016962    6
0   101_DEN084Y5_15_E05_010_clean.fcs   (root,) Time    RectangleGate   None    None    284846  99.844369   99.844369   1
1   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time)    Singlets    PolygonGate None    Time    235331  82.488345   82.616923   2
2   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets)  aAmine- PolygonGate None    Singlets    160955  56.418031   68.395154   3
3   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-) CD3+    PolygonGate None    aAmine- 131305  46.025097   81.578702   4
4   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+)   CD4+    PolygonGate None    CD3+    81135   28.439483   61.791249   5
5   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+)   CD8+    PolygonGate None    CD3+    46546   16.315328   35.448764   5
6   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) CD107a+ RectangleGate   None    CD4+    70  0.024536    0.086276    6
10  101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) CD107a+ RectangleGate   None    CD8+    306 0.107259    0.657414    6
7   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) IFNg+   RectangleGate   None    CD4+    9   0.003155    0.011093    6
11  101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) IFNg+   RectangleGate   None    CD8+    348 0.121981    0.747647    6
8   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) IL2+    RectangleGate   None    CD4+    8   0.002804    0.009860    6
12  101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) IL2+    RectangleGate   None    CD8+    65  0.022784    0.139647    6
9   101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD4+) TNFa+   RectangleGate   None    CD4+    28  0.009815    0.034510    6
13  101_DEN084Y5_15_E05_010_clean.fcs   (root, Time, Singlets, aAmine-, CD3+, CD8+) TNFa+   RectangleGate   None    CD8+    336 0.117775    0.721867    6


text:
How does one distinguish nuclear DNA from mitochondrial DNA when doing WGS?
I'm interested in doing de-novo sequencing but also phylogenetic analysis. In particular, after de-novo sequencing and annotating the genome, I need to align the CO1 gene and the nuclear 28S rRNA gene of several species. When extracting DNA and sequencing with PacBio for example, how does the assembler know what corresponds to nuclear and what to mitochondrial DNA?

Does it need to? I mean, and I stress that I have never done de-novo sequencing so I may be completely wrong here, but won't the assembler simply try to assemble the longest contiguous contigs it can? And, if so, won't that naturally result in the mitochondrial sequences getting assembled together in one set of contigs and the non-mito sequences in others? I would expect this to be even easier with long read technology like PacBio. Finally, if you're only interested in a couple of genes, can't you directly target the highly conserved parts of them to amplify and sequence?

Cannot comment on CO1, but rrRNA should be different enough between the mitochondria and the nucleus as to not be confused. For one, mitochondria (like bacteria) should have 23S subunits instead of 28S. Many programs (e.g. barrnap, whatever is used by Prokka) automatically distinguish them.

Could you kindly review the answer provided? Please do consider upvoting or accepting an answer.

According to this pdf,  from pacb.com
isolation with mitoISO MITOISO
The kit will include three steps: 1) mechanical
rupturing of tissues or cells 2) removing cellular debris and nuclei by low speed centrifugation and 3) harvesting
mitochondria by high speed centrifugation.
Then, you can follow the DNA extraction with mitochondrial DNA, and then sequenced on the PacBio Sequel.
If you are doing whole genome sequencing, the mitogenome will have a higher read depth according to this article stating that lower read depth of the nuclear genome separating the two different origins:

those of mitochondrial origin & those of nuclear origin


text:
Is there an alternative to bulked segregant analysis for insects?
This strategy seems to be most commonly used for plants. When crossing animals isn't possible, how can I do a similar study to identify a particular locus responsible for a polymorphic trait (specifically, the pigmentation of a beetle)?

I'll have a go at answering this. I don't understand beetles within Insecta: I am well aware they very distinct, obviously in phenotype and ecology, but also genetically. I think they have unique compositional bias even amongst holometabolous insects.
Personally I think if you've a recent genetic trait, e.g. insecticide resistance, that evolved within the last 30 years or less, you're dealing with pretty homogeneous backgrounds between traits. In this example trait, bulk segregant analysis might work.
Beetles are really ancient. I suspect their families are particularly ancient, because they are just so successful and I'd guess Scarabaeidae are ancient too, i.e. the family scarab beetles. This family could easily be 50 million years old, or older. Each species therein could a few million years old. For insects thats normal and these are beetles.
The concern here is over million+ years of evolution (I don't know the study size), they will comprise loads of neutral signal that would act as noise, i.e. there will be excessive third codon wobble which will mask SNPs/genes causing colour (color) variation.
My approach, given I know nothing about colour (color) variation in insects other than "white-eyed" phenotype in Drosophila melanogaster,

assess genetic complexity of colour variation in other insect Orders;
get hold of as many proteins involved in colour in insects as possible;
particularly get hold of colour variation proteins in beetles;
sequence the genome of a representative scarab beetle (which I think your doing)
perform psi-blast using the hypotheses from 1., 2. and particularly 3. to assess which proteins in the scarab beetle might be responsible.
perform deep sequencing on those targets, because they could be multigene families *.

The next step is phylogenetics to assess convergent and/or parallel evolution **. That in my view is a separate question because its quite complicated and you'd need to write your own pipeline to automate it (you could do it by hand however).
What I would be cautious about is that beetles are not any other insect so inter-Order homology needs to be treated with caution.

*, Just to point out, there are lots sequencing strategies its not my thing, e.g. low coverage draft genomes. PCR deep sequencing for a given gene target, without any draft genomes and isn't conserved e.g. rRNA types genes, will need small PCR fragments because of the neutral mutations which will disrupt primer specificity. If you're good at sequencing a draft genome could be cool (just guessing).
**, A phylogenetics approach would dovetail with the COI and rRNA gene sequencing (28S?) already underway, because these genes should approximate to the species phylogeny. What is sought are residues  within the colour genes (or perhaps presence/absence genes) that will cluster not according to their evolutionary descent but according to their colours (colors).


text:
Fixing FASTA file for Local BLAST Database
I recently prefetched 157 SRA files from an NCBI BioProject using the SRAtoolkit. I then used the toolkit to download those files in FASTA format.
Each individual FASTA file looks something like this:
(example for SRR5678544 file)
>SRR5678544.1 HWI-ST1146_88:5:1112:6472:81473 length=194
CCATGCGGGGTATCGTATGCTTCCTTCTGCACTACCCTTTAGCTGTTCTATATGCTGCCACTCCTCAATTGGATTAGTCTCATCCTTCAATGCTATCACAAGAAGTAGAGACACAAATGCAAGAGGAGCATATAAATTACAAAACACCATCACTGAGGGCCCTAAAGCGGTTCCCACGAAAAAAAGGAGAGTAG
>SRR5678544.2 HWI-ST1146_88:5:1113:13218:62635 length=194
 TGGATGGTGGTTGGGAGTGGTAAGGTTGAATGAGACACGGTAACGAGTGGGGAGGTAGGGTAATGGAGGGTAAGTTGAGAGACAGGTTCGTCAGGGGACACACCACACCACACCCACACCCACACACCCACACACCACACCCACACACCACACCCACACCCACACCCCCACACACCACACCCCCACCCCACCCA
>SRR5678544.3 HWI-ST1146_88:5:2206:11224:22269 length=194
 GTCTCTTAACTTACCCTCCATTACCCTACCTCCCGACTCGTTACCCTGTCTCATTCAACCATACCACTCCCGACCACCATCCATCCCTCTACTTACTACAGTATGGTGAGTGGGACATGGTGGATGGTAGGGTAAGTGGCAGTGGAGTTGGATATGGGCAATTGGAGGGTAACGGTTATGGTGGACGGGGGGTG

And so on and so forth. (Ignore the " and ' marks I had to input those to format it right on here)
Question 1: Can someone explain why this file has so many sequences (>)? Every FASTA file I have worked with has one line with ">" at the top and a header sentence followed by a single big paragraph of the genome sequence.
Question 2: If I wanted to combine all the sequences in a single FASTA file to look like a "normal" FASTA file - just a single paragraph - how would I do that?
Reasoning: I'm asking simply because I need to make a BLAST database out of the 157 sequences and the "makeblastdb" command only takes in a single file as the input. I was going to try and combine each individual FASTA file into single paragraphs and then make a huge FASTA file that includes all 157 sequences separated by ">".

This is a fairly straight forward question, generally, someone here should answer this, but your biological objectives are quite important. Firstly, do you have a reference genome, or is that going to be the purpose of the Blast? What are you looking for in the Blast, is it some sort of mapping?

I've edited the sequences into a coding block, rather than raw text. When I did this the taxa id (Genbank code and data) were on the same line as the sequence data. There was a carriage return this line and the next one starting with  >. Please check whether my editing is what was in the original file, i.e. I inserted a carriage return between the taxa id, e.g. >taxa id and the sequence, e.g. GTCTCTTA. Thats quite important

crossposted https://www.biostars.org/p/9537532/

The question appears to be a confusion regarding short read data. The point fastqc is normal for assemblies rather than raw fasta is a good point, ie. its not cool to make an assembly directly from the data set in the question.

I was confused about the data but now I understand it's short unassembled reads. My next question is then what now? Do I assemble them? Or do I align them? I have a reference genome I could use but the purpose of the BLAST database is to eventually do some phylogenetic analysis, analyzing closely related species, looking at specific genes in those species, etc.. There's a lot I want to be able to do with the data I am collecting I am just confused on how to process them before those analysis.

For instance, creating a local BLAST database requires a SINGLE fastA file as the input. And as I have 157 files (to start) I need to somehow combine all the 157 sequences into one file. I can use this line of code to acheive that: cat *.fasta > allfiles.fasta but I can't include all of the short raw reads without making them a consensus or aligned or assembled sequence first. That's the part I'm stuck at mainly.

Briefly the best strategy I think is to download fastQ files and process this via an assembler. I wouldn't enter raw short read Illumina data into a local Blast database. What I'm pretty certain of is thats a bad move.
Which assembler I personally think is a separate question on Bioinformatics SE. It will depend on the genome, for large eukaryote genomes BWA-mem2 seems popular, but its not my thing and there's other stuff upstream and downstream (someone might dispute that). For assembling small genomes, Spades is popular, but preferences are changing and there's also upstream QC. The other question is whether you have a reference genome to assist assembly.

The update is this a yeast genome with a reference. This is complicated because its a small genome but will comprise heterozygosity. If it was me - i.e. someone way outside the yeast community  - I'd use Spades and provide a reference genome to perform an assembly. However, I'm not sure thats at all cool, because you will will reduce the eukaryotic genome to a haplotype. This could get complicated because what you would be left with may not be biologically real, e.g. two alleles at different parts of the same chromosome on different strands may result in a in silico artefact. I think its a separate question personally.

I do have a reference genome (S288C: https://www.ncbi.nlm.nih.gov/nuccore/NC_001224.1?report=fasta) that I was planning on using. However, I would like to implement de novo assemblies as well just to have that in case it's needed. I have some experience with SPAdes but I assume the tool used doesn't really matter too much?

My other question is why would you assemble over aligning? These are all yeast genomes so they are pretty ambiguous in terms of their ploidy/heterozygosity/etc., so I want to make sure I'm doing the right thing for my data. This question is purely curiosity driven and I really appreciate your help so far.

Hi @rimo on Bioinfo SE this is a separate question. For me its a difficult question because its a eukaryote but a small genome. The issue is (I presume) there'll be heterozygosity in the data and I don't know how the yeast community deal with this (they might ignore it). Most small genomes (my stuff) are haplotypes ( obviously I would ignore it). Most higher eukaryotes use VCF, but this is yeast. It gets complicated because you can't work out the linkage so decomposing a eukaryote into two haplotypes isn't possible. Personally I would assemble via a reference genome. I've edited the reply.

That's exactly what I was thinking. I just posted a separate question to the Bioinfo SE about this exact thing. Thank you so much for all your help with everything else. I'm going to try and assemble via a reference genome to start and then go from there.


text:
Advise on building an effect ML model for predicting important proteins for drug response
I want to create a model to predict proteins which could be associated with drug response in cancer cell lines. I have cell line proteomics data, with compound screening data they have gained for a number of these cells. My aim is to relate these proteins to drug response. I am using this exercise as an introduction in the use of ML in biology.
To this end, I have created dictionaries, with selected compounds as the keys, and their respective datasets. Cells are rows, proteins are columns, and the drug response is the y variable.
I have approached the task as a regression task. The aim being to predict the ln(IC50) - a measure of the concentration required to induce cell death - and subsequently pull out the features (proteins) which have made the largest contribution to this prediction. As my baseline model, I have attempted linear regression. However, my model performance has consistently been very poor. A problem with this dataset is the large number of missing values (heatmap below - yellow = missing). However, even when trying a simple feature selection strategy with sklearn (first function), the explained variance is around 0.3 for the best performing compound. I have plotted the residuals in attempt to dig into this further, but the spread looks quite distributed.
Any advise/suggestions would be appreciated. Or any further reading/papers relating to ML in biology which could help a beginner build and use effective models.
Training and test functions shown below.
    def process_training_and_test_data_for_LG(df):
        training_datasets = {}
        test_datasets = {}
        for drug in df.keys():
            df[drug].dropna(axis=1,inplace=True)
            print("at "+drug)
            X = df[drug].drop(["cell_line_name", "ln_IC50","putative_target"], axis=1)
            y = df[drug]["ln_IC50"]
            X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)

            feature_selector = SelectKBest(score_func=f_regression, k=50)
            feature_selector.fit(X_train, y_train)

            X_train = feature_selector.transform(X_train)
            X_test = feature_selector.transform(X_test)

            #Add data to relevant dictionaries
            training_datasets[drug] = [X_train, y_train]
            test_datasets[drug] = [X_test, y_test]
        return training_datasets, test_datasets

def train_models(training_data):
    trained_models = {}
    for drug in training_data.keys():
        #initiate model
        lm = LinearRegression()
        #extract drug specific training data
        X_train = training_data[drug][0]
        y_train = training_data[drug][1]
        #fit model
        model = lm.fit(X_train, y_train)
        #Append dict of trained models
        trained_models[drug] = model
        print("Model for {} trained".format(drug))
    return trained_models

def test_models(model_dict, test_data):
    prediction_dict = {}
    for drug in model_dict.keys():
        drug_model = model_dict[drug]
        X_test = test_data[drug][0]
        predictions = drug_model.predict(X_test)
        prediction_dict[drug] = predictions
    return prediction_dict

The missing values might be an issue indeed. You might want to use imputation methods, e.g.
from sklearn.impute import KNNImputer
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer

def impute_df(df, feature_columns_to_apply, strategy):
    if strategy == 'mean':
        for col in feature_columns_to_apply:
            df.loc[:, col] = df.loc[:, col].fillna(df[col].mean())
    elif strategy == 'median':
        for col in feature_columns_to_apply:
            df.loc[:, col] = df.loc[:, col].fillna(df[col].median())
    elif strategy == 'knn':
        knn_imputer = KNNImputer(n_neighbors=5)
        df_tmp = pd.DataFrame(knn_imputer.fit_transform(df[feature_columns_to_apply]), columns=feature_columns_to_apply)
        df.loc[:, feature_columns_to_apply] = df_tmp[feature_columns_to_apply]
    elif strategy == 'iterative':
        iterative_imputer = IterativeImputer()
        df_tmp = pd.DataFrame(iterative_imputer.fit_transform(df[feature_columns_to_apply]), columns=feature_columns_to_apply)
        df.loc[:, feature_columns_to_apply] = df_tmp[feature_columns_to_apply]
    return df

df_all = impute_df(df_all, feature_columns_to_apply, 'mean')

(IterativeImputer is heavy and slow)
These are just sample strategies, you might want to use other methods, see e.g. here.
Also consider other models such as random forest / XGBoost, they outperform linear regression in many cases.

The gist looks fine ..

pandas - looks fine
SelectKBest - cool
transformation (obviously there's a lot of them)
standard ML

A common ML situation is, the classifier working but the accuracy is notably lower than expected. My personal first thing to check is the use of linear regression as a baseline model. The scatterplots give a gist, we don't know in particular what data points belong to what classification - thats crucial. What's helpful are coloured dots to assess how they sit against the regression line.
What I think has happened is there's been a PCA via SelectKBest, unsupervised learning, which doesn't account for much of variance (cf 0.3), i.e. this is not an accuracy score.
I've three suggestions:
1/ An acceptable baseline model often used is k-nearest neighbour ..
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
clf = KNeighborsClassifier(n_neighbors=3)
clf.fit(X_train, y_train)
print("Test set accuracy: {:.2f}".format(clf.score(X_test, y_test)))

I think this will produce much better results with an accuracy of >0.7
2/ However,  I would check the following (please see caveat),
for drug in df.keys():
    df ....
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)

Just make sure that X_train, etc ... is supposed within the loop, rather than immediately outside it (four spaces to the left). Looking at your outputs they look fine.
The thing to check is its not:
for drug in df.keys():
    df....
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)

The above loop makes sense to me.
The original loop position train_test_split looks like the key loop might repeatedly over-ride the X_train, X_test, y_train, y_test outputs.
3/ Very final point SelectKBest followed by fit. What about fit_transform instead? If the issue is the features are resulting in low description of the variation, maybe thats key here?

Caveat I do stress that its hard to comment beyond general terms because a good biological understanding is needed of the goal of the ML training. In context its not my area and I haven't read any of the  papers.

Many thanks for this. Apologies there was an error in the indentention in the code I posted, I have amended this in the edit.

Regarding the the use of KNN, I should of mentioned I have approached this as a regression task, and the y variable is a measure of the concentration required to induce cell death) I'll update this. Though changing this to a classification task could simplify it and allow me to look at feature contribution to classification.
For the fit/fit_transform point, I thought that fit_transform did the same function performed the same as .fit() followed by .transform(), is there a subtle difference?


text:
Single or partitioned model for supermatrix tree inference?
I am dealing with several genes in my dataset. For each of the genes, I have built gene tree, with their best estimated model. I am intended to look into the effect of the tree from concatenated sequences (supermatrix tree), using Maximum Likelihood approach.
I am considering two options to apply the model for the supermatrix tree inference.

Obtain a single best estimated model for the supermatrix gene
Apply partitioned model for the supermatrix gene, based on each gene's model

From what I have read in most if not all of the research articles, they seemed to be using option 1, where a single model applied for the supermatrix tree inference. Since the purpose of using evolutionary model is to describe the changes in the gene and given that the evolution pattern in each gene differ, why wouldn't the respective model of each gene be considered in the tree inference?
For examples,
Dataset A: gene1 and gene2 are having LG and WAG model respectively as the best estimated model. But when the two genes are being concatenated, it turned out that LG is the best model for the matrix. It seemed that most if not all of the experiment opts for option 1 in this case. Theoretically, will the supermatrix tree using option 1 and option 2 in the model application still similar? Also, will the evolutionary pattern in gene2 still being represented correctly on the tree when we opt for option 1 in the supermatrix? Or it actually does not impact much on the calculation of the tree?
Dataset B: gene1 and gene2 are having LG as their best estimated model. The concatenated matrix has LG as its best estimated model too. Will model partitioning has any effect in this case? Does the ML tree/site likelihood calculation affected by the length of the input sequences?

I understand concatenation well. What I'd like to be clear about is the 'super-matrix', is this the single mutation matrix derived from the concatenated alignment? Secondly are you obtaining your matrix before you perform the tree calculation and then inputing this matrix for the ML tree search?

@M__ yes for both. Say if I have 5 genes. The supermatrix derived from adding gene 2 at the end of gene 1, gene 3 at the end of gene 2, ... , gene 5 at the end of gene 4. They will exist as a big and single concatenated alignment matrix, concatenating gene1 to gene5 side by side. Then the supermatrix would be used together with the model/s to calculate for the ML tree.

Thanks, two final questions are these genes physically linked? Or random amplicons, across a genome? Finally, is this a mix of rRNA genes and protein genes and within the protein genes are they all single copy house keeping genes? I'll answer before the end of the week.

The genes are random across genome and yes they are  all single copy rRNAs and protein genes.

Thanks. Sorry to be persistent, are the protein genes 'house-keeping', e.g. often a core metabolic enzyme and all part of the core-genome, specific to the Order/family/species or a mix of the two?

No worries. It's just that my question is actually referring to those in general situation. It would be great if you could help me understand the concern. Probably I could justify it clearer.

Okay, you can't concatenate genes into a singular matrix particularly if these are rRNA genes versus protein genes, except in rare situations where the tests say 'thats okay' (below). This is simply because they evolve at very different rates rRNA genes have loci which evolve very slowly and are highly conserved.
The theory goes back to incongruence length difference test (ILD), which is somewhat confusingly also known as the partition homogeneity test and was then modified to become the localised ILD test.
The modern version of this test, simply called 'homogeneity' is found in IQ-Tree. Thus, the best approach is to test the assumption empirically,

IQ-TREE provides three matched-pairs tests of symmetry (Naser-Khdour
et al., 2019) to test the two assumptions of stationarity and
homogeneity. A simple analysis:
iqtree2 -s example.phy -p example.nex --symtest-only

Just to be clear this tests whether the rate between two partitions (genes) is homogeneous. If 'yes' the supermatrix is fine, if 'no' the data needs partitioning.
The theory is simply that if two genes evolve at noticeably different rates they can't be combined into a singular matrix. The theory goes back a long way is used to circumvent long-branch attraction issues. Thus, at a minimum the rRNA genes and protein genes need to be compartmentalised, i.e. have their own matrix. Between protein genes, if they are house-keeping genes, generally they'll have homogeneous rates (which you can empirically test via IQ-TREE).
If someone is concatenating rRNA with protein genes without accommodating rate variation via partitioning it is simply sloppy phylogenetics and if published represents a shortcoming with the reviewing process. This is because in all certainly these loci will fail any rate homogeneity test, there's a slim off-chance its okay which you can test it empirically. Say if 18S and 28S were concatenated (or just 16S if these are bacteria) with an atp gene, that can give a different tree to 18S with atpA, atpB, atpC, atpD etc... so the ratio between rRNA to protein genes can start influencing the tree ... which is simply wrong.
If someone concatenates a load of house-keeping genes thats okay because there's a good chance they're homogeneous and if even the test failed, they may not be that far out, rRNA versus proteins can be way out. This can easily be an order of magnitude out and thats difficult to ignore. When there's proteins under selection versus conserved housekeeping genes ... its very unlikely homogeneity would hold.
Summary As a general rule of thumb in your method of tree building compartmentalising all genes is safer route, computationally there's no difference between zero partitions and a partition for every gene. You are never going to make a mistake.
Where it gets complicated is when the phylogenetic parameters, particularly in advanced models, are themselves estimated via maximum likelihood alongside the topology. The rational for this is complicated. Computationally, this can get very, very expensive even on supercomputers and hence there's a tendency to make an assumption (often without checking it) that there are fewer partitions than there actually are. Regardless, assuming away rRNA versus protein genes even in this method of tree estimation should not done and the test will almost certainly fail.

Appreciate your answer. I agree with the possible induced error when concatenating together rRNAs and protein genes. Sorry for the confusion caused from the comments. What I mean was I have dealt with both type of datasets and they were concatenated separately. But that was quite an interesting input for me, especially the ILD test, will try that out.

Btw, from your last paragraph, though not really familiar with the concept behind, I was aware that it is possible to estimate best models via the ML algorithm. But what does it meant by "there's tendency to make an assumption that there are fewer partitions than there actually are"? Also, I am very interested with "computationally there's no difference between zero partitions and a partition for every gene", would you mind to elaborate a bit more if possible? What do you think about it theoretically? Btw, I have updated the question. Thanks!

Agreed @web the error was unfortunate. "computationally no difference ...", an exaggeration, should be "the difference in computation time is not a major issue". Final point "tendency ...", a reviewer will be aware of the difference between what is computationally practical against the theoretical requirements of the calculation, certain approximations are considered ok.


text:
snRNA-seq Demultiplexing before or after QC
I'm currently analyzing snRNA-seq data and have to demultiplex them using HTODemux in Seurat.
The question is : should I do my QC filters on my droplets before or after the demultiplexing ?
From what I can test, the QC filters actually is going to have an effect on the assignation to each sample, as you can see in the small summary data below
Demultiplex, then QC :
Cutoff for 2-2 : 104 reads
Cutoff for 2-3 : 242 reads
Cutoff for 5-2 : 80 reads
Cutoff for 5-5 : 78 reads

 Doublet Negative      2-2      2-3      5-5      5-2
     347     1754     1182      960     1278      959

QC, then Demultiplex :
Cutoff for 2-2 : 109 reads
Cutoff for 2-3 : 194 reads
Cutoff for 5-2 : 83 reads
Cutoff for 5-5 : 74 reads

 Doublet Negative      2-2      2-3      5-5      5-2
     426     1466     1066     1284     1302      936

As you can see, in the case of sample 2-3 we have a 25% difference in assignation, so I don't think we can completely ignore that.
To be clear, I understand why there's a difference (normalisation, clusterisation differences, quantiles differences), I am looking for best practices here.
Thanks in advance

In the end I took the decision to first do the QC and then Demultiplex, reasonning that the Demultiplex steps (normalisations, clustering etc ...) are probably cleaner this way, leading to a (I hope) better assignation


text:
correlation between imputed genotype and true genotype
There is a graph in the link https://odelaneau.github.io/GLIMPSE/rsquare_eur.html. This graph shows the imputation accuracy of certain MAF bins. Now I understand that $r^2$ is the correlation between imputed genotype and true genotype. What I don't understand is how they calculate this $r^2$ correlation. Can anyone suggest to me how I can define the r^2 correlation between true genotype and imputed genotype? I suppose it depends on the individual who defines his own $r^2$. So, it may differ for the different algorithms. However, different aspects or different $r^2$ definitions will help me understand and represent this graph slightly differently.

This is done by downsampling. Take the 1000 genomes, set some genotypes as missing ./., impute them using GLIMPSE, then measure correlation between the genotype dosages in the imputed dataset and the same sample at high coverage - in this case, the 1000 genomes sequenced to 30x coverage by the NY Genome Centre.

Is the correlation a Pearson correlation coefficient?

@winni2k judging by this yes I think so https://github.com/odelaneau/GLIMPSE/blob/master/concordance/src/utils/basic_stats.h#L110


text:
What is the best way to process yeast genomes?
I have obtained several hundred raw, unassembled yeast genomes from NCBI and I am looking for advice on how to process the genomes for downstream analysis.
I have a reference genome (S288C) to use for assembly/alignment but I'm not sure what to do first.
So far I have run a quality assessment (FastQC) followed by trimming (Trimmomatic). But I'm not sure if I should align or assemble next. The problem with yeast is if I assemble to a reference genome it reduces the genome to a haplotype and could miss moments of heterozygosity present. And I'm not sure how it would work to do a de novo assembly as I'm assuming it's the same problem. Same with aligning to the reference genome.
Are there ways to work around this? Does the yeast research community just ignore this issue? Are there tools currently out there that take all of this into account?

What is your data type? Many assemblers can resolve heterozygosity.
PacBio HiFi reads can be used with hifiasm to resolve both haplotypes rather trivially. This is the best case scenario and older sequencing projects are unlikely to have HiFi reads.
If they are merely PacBio subreads, you can use Falcon-Unzip to get heterozygous regions (though this is now quite old and I am not sure that it is recommended anymore).
For Illumina, the options are worse. Platanus is an assembler intended for the purpose of assembling heterozygous genomes. If nothing else, you can make the assembly and then just call variants.
Note that many yeast are probably sequenced from a haploid strain.

I have both PacBio HiFi reads as well as Illumina seqeunces. I'll defineityl try out hifiasm for the HiFi reads but I have a ton of downloaded Illumina sequences that I don't know how to process.
If I were to perform a de novo assembly with SPAdes on those sequences and then called the variants with bcftools or something of the sort would that give me reliable data?

@rimo looking in more detail, there are better Illumina assemblers for highly heterozygous genomes. Updating answer with this info.


text:
Naming conventions for kegg nodes with multiple genes
So I'm downloading XML files from Kegg using BioConductor and and I'm running into a problem. If we look at HIF-1 signaling pathway we see the node growth factor (GF) on the left hand side, we can trace it's XML entry to this entry:
    <entry id="17" name="hsa:1950 hsa:3479 hsa:3630" type="gene"
        link="https://www.kegg.jp/dbget-bin/www_bget?hsa:1950+hsa:3479+hsa:3630">
        <graphics name="EGF, HOMG4, URG..." fgcolor="#000000" bgcolor="#BFFFBF"
             type="rectangle" x="87" y="272" width="46" height="17"/>
    </entry>

While there is no GF on this, there are the three info card links that I can use to extract information. However, when I click on that link, there isn't a place on those info cards that says that these three genes are related to that general title of GF.
Normally this isn't a problem, I would just look at the map, but I'm trying to do something at scale where I can't look at every instance of this.
So my question is, how do I find the labels that Kegg uses in its map using the XML file/python package?

Python's Bio.KEGG.Gene has an attribute called name described as "A list of the gene names". Might help? Otherwise "orthology" is its own universe. Last ditch attempt try beautifulsoup and webscraping. In theory webscraping will work because if you can access the links via clicking them beautifulsoup will do this too.

Thanks, I'll check out the name attribute! I'm trying to avoid scrapping at all costs cause this particular use case won't be pretty, but if it's all there is then...

How do you download the files from KEGG with Bioconductor? It might be easier to use Bioconductor to retrieve the gene names

If you used KEGGREST you should be able to retrieve the gene names of each pathway: https://bioconductor.org/packages/release/bioc/vignettes/KEGGREST/inst/doc/KEGGREST-vignette.html

The problem with the name attribute is that it returns multiple names for each card inside that node. GF for example returns insulin, insulin like growth factor 1, and epidermal growth factor. Which one do I use? and how did the original map just label it as GF (growth factor)?


text:
What data set was used Gene Expression Data of Postmortem Tissues?
I want to run the experiments "Sample Application to Genotype-Tissue Expression (GTEx) Project Gene Expression Data of Postmortem Tissues" mentioned in "Enter the Matrix: Factorization Uncovers Knowledge from Omics" by Stein-O'Brien et.al link
The paper metioned

dataset from postmortem samples in the GTEx project [42] (Figure 3).
Specifically, we select a subset of GTEx data containing 12 brain
tissues in the GTEx data for seven individuals for which we had
previously performed only NMF analysis [41] (codes are provided in the
supplemental information online).

I check with https://gtexportal.org/home/datasets and assumed that the paper used the GTEx Analysis V6.
But which files do I need to download? I am computer scientist I am unsure which ones contain the brain  tissues
Further I found
Brain_Anterior_cingulate_cortex_BA24_Analysis_cis-eQTLs.txt.gz  1.5G
Brain_Caudate_basal_ganglia_Analysis_cis-eQTLs.txt.gz   1.9G
Brain_Cerebellar_Hemisphere_Analysis_cis-eQTLs.txt.gz   1.7G
Brain_Cerebellum_Analysis_cis-eQTLs.txt.gz  1.9G
Brain_Cortex_Analysis_cis-eQTLs.txt.gz  1.9G

But I am unsure which one was used in the paper


text:
Mass spectrometry: How to obtain expression counts for protein
I was given a data set consisting of 6 samples: 3 control and 3 treated. I dont have information about how data were generated and what the outputs represent.
Each sample corresponds to a folder within which there are the following files:
"DB search psm.csv" "de novo only peptides.csv" "de novo only peptides.xml" "de novo only spectra.mgf" "filter.csv" "peptide features.csv"
"peptide.csv" "peptides_1_1_0.mzid" "peptides.pep.xml" "protein-peptides.csv" "proteins.csv" "proteins.fasta" "PSM ions.csv"
The "peptide.csv" files have the following column names:
Peptide -10lgP Mass Length ppm m/z RT 1/k0 Range Area Ctrl-2 Fraction Precursor Id Source
The "proteins.csv" files have the following column names:
Protein Group Protein ID Accession -10lgP Coverage (%) Coverage (%) Ctrl-1 Area Ctrl-1 #Peptides #Unique #Spec Ctrl-1 PTM Avg. Mass Description
Has anybody ever worked with these type of data?
I would like to do exploit these data and perform a differential analysis but i don't know how to create a SummarizedExperiment from these data.
Any help would be greatly appreciated.
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.

Could you provide an example dataset? And post the headers and top 5 rows of each file?


text:
I want to map TFBS on to the sequences from dataset containing 8448 (210 base) sequences. its not giving all the hits I need
I want to map TFBS (transcription factor binding sites) on to the sequences from dataset containing 8448 (210 base) sequences, but it is  not giving all the hits I need
import pickle
from pickle import load, dump
from collections import defaultdict
from Bio.Seq import Seq

#first element is BREu = SSRCGCC, length = 7 bases, size of matrix = 4*7, max score = 3.33
# first row = A, second row = T, third row = G and forth row = C
BREu = [
    [0.0, 0.0, 0.5, 0, 0, 0, 0],
    [0.0, 0.0, 0.0, 0, 0, 0, 0],
    [0.5, 0.5, 0.5, 0, 1, 0, 0],
    [0.5, 0.5, 0.0, 1, 0, 1, 1]]

#second element is ETS = MGGAART, length = 7 bases, size of matrix = 4*7, max score = 6
ETS = [
    [0.5, 0, 0, 1, 1, 0.5, 0],
    [0.0, 0, 0, 0, 0, 0.0, 1],
    [0.0, 1, 1, 0, 0, 0.5, 0],
    [0.5, 0, 0, 0, 0, 0.0, 0]]

#Third element is INR = YYANWYY, length = 7 bases, size of matrix = 4*7, max score = 3.75
INR = [
    [0.0, 0.0, 1, 0.25, 0.5, 0.0, 0.0],
    [0.5, 0.5, 0, 0.25, 0.5, 0.5, 0.5],
    [0.0, 0.0, 0, 0.25, 0.0, 0.0, 0.0],
    [0.5, 0.5, 0, 0.25, 0.0, 0.5, 0.5]]

#forth element is NFY = CCAAT, length = 5 bases, size of matrix = 4*5, max score = 5
NFY = [
    [0, 0, 1, 1, 0],
    [0, 0, 0, 0, 1],
    [0, 0, 0, 0, 0],
    [1, 1, 0, 0, 0]]

#seventh element is SP1_NC = GGGNGGG, length = 7 bases, size of matrix = 4*7, max score = 7
#Gives some possibilities of having any bases at position 4 while still favoring G
#threshold score >= 6.5 : keep all degenerated motif with degenerated = 6.5 and canonical = 7

SP1 = [
    [0, 0, 0, .5, 0, 0, 0],
    [0, 0, 0, .5, 0, 0, 0],
    [1, 1, 1, .5, 1, 1, 1],
    [0, 0, 0, 1., 0, 0, 0]]

#sixth element is TATA = TATAWAAG, length = 8 bases, size of matrix = 4*8, max score = 7.5
#threshold score > 20 : keep all degenerated motif with degenerated = 20 and canonical = 23
TATA = [
    [0, 4, 0, 4, 4, 1, 1, 0],
    [4, 0, 4, 0, 4, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 1],
    [0, 0, 0, 0, 0, 0, 0, 0]]

# define matrix and thresholds as tuple
motifs = dict(
    BREu = (BREu, 5.5),
    ETS = (ETS, 6),
    INR = (INR, 3.75),
    NFY = (NFY, 5),
    SP1 = (SP1, 6.25),
    TATA = (TATA, 20))

def get_motif_score(seq, motif):

    motif_size = len(motif[0])
    bases_encodings = {'A': 0, 'T': 1, 'G': 2, 'C': 3}
    coded_seq = [
        bases_encodings[b] if b in 'ATGC' else None
        for b in seq]

    return [
        sum(
            0 if coded_b is None else motif[coded_b][i]
            for i, coded_b in enumerate(coded_seq[start:start+motif_size])
        ) for start in range(len(seq) - motif_size)]

def get_motifs(seq):
    """find motif occurences in given sequence

    return list of (motif sequence, start, end, score)

    """
    seq_size = len(seq)
    motif_hits = {}

    for motif_name, (motif_matrix, threshold) in motifs.items():
        motif_size = len(motif_matrix[0])
        rcseq = str(Seq(seq).reverse_complement())
        scores = get_motif_score(seq, motif_matrix)
        rc_scores = get_motif_score(rcseq, motif_matrix)

        motif_hits[motif_name] = [
            (seq[i:i+motif_size], i, i+motif_size, score,'+')
            for i, score in enumerate(scores)
            if score >= threshold
        ] + [
            (rcseq[i:i+motif_size], seq_size-(i+motif_size), (seq_size-i), score,'-')
            for i, score in enumerate(rc_scores)
            if score >= threshold
        ]

    return motif_hits

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.

What is the use and how we should interpret the matrix in your code ???

Hi pippo1930, The matrices are self designed to map the motifs on basis of arbitary score, with the function of getting score, I am calculating the score on the window size equal to the length of the motif (length of matrix row), then if the score is above the threshold as indicated, it means that the motif exist. note that the sequence are non - consensus, which means that some basis are not very particular (that is why the score to map that particular base is less). I am not sure if that answer your question and thank you for taking interest to solve this problem.

Its a really hard alghiritm to understand, here usually you should provide a test input with the expected results vs actual ones to help people understand the issue

How of curiosity could you show me how a BREu motif can get a > 5,5 score

print(get_motifs('TATAWAAG')) gives {'BREu': [], 'ETS': [], 'INR': [], 'NFY': [], 'SP1': [], 'TATA': []}

lower the threshold ???

or change the matrix values


text:
Where can I download 30x 1000 genomes cram files?
From the preprint published by 1000 genome project (https://www.biorxiv.org/content/10.1101/2021.02.06.430068v1.full) I think the 30x data is for WGS. Can anyone confirm for me if the following file is really WGS data?
ftp://ftp.sra.ebi.ac.uk/vol1/run/ERR324/ERR3242643/HG03007.final.cram.
Also, Is there any more freely available sequencing data where I can find WGS data for 30x and low coverage(<10x) for the same sample for the South Asian population?

Yes, those are the .cram files for the 30x whole genome sequencing data for 1000 genomes.
If you want 10x and 30x sequence data for the same individual, you can just downsample the 30x .cram file using samtools view -s.

30x WGS or WES? Can u specify it in your answer?

WGS - I don't know why you'd think it was exome sequencing - there's no mention of that in the paper or the website.

There is actually. But that was at the depth of 65.7x coverage. However, I asked just to be sure, anyway once again thanks for your answer.


text:
Creating multiple phenotype datasets using bootstrap method "Bootstrap-samples-by-column-of-a-data-frame-in-r" for DEG analysis
I am working on a datasets and after some discussion with my group, we doubt that maybe one or more of our controls are different than the other controls. The motivation is to see if one or more controls have been effected differently by the solvent they were kept in.
I have been suggested to use bootstrap method. If we suppose that I have a dataset with 5 controls and 3 treated samples. I want to create 5 new data frames with the information such as every new data frame skips one out of 5 controls and do resampling with replacement. I want to see how stable the controls are and by skipping a specific control, would the DEGs result change.
Let us suppose that the original data frame is like this:
x <- round(matrix(rexp(480 * 10, rate=.1), ncol=8), 0)

rownames(x) <- paste("gene", 1:nrow(x))

colnames(x) <-c("control1","control2","control3","control4","control5","treated","treated","treated")

I want to create 5 new data frames (as there are 5 controls in this study) where each data frame skips one specific control and replace with with another control (which means some other control will repeat).
For example one of the 5 data frame can look like:
x1 <- round(matrix(rexp(480 * 10, rate=.1), ncol=8), 0)

rownames(x1) <- paste("gene", 1:nrow(x1))

colnames(x1) <-c("control1","control1.1","control3","control4","control5","treated","treated","treated")

You can see that this new data frame skipped control2 with a copy of control1 called control1.1.

we doubt that maybe one or more of our controls are different than
the other controls.

Are they distinct by, say, PCA?


text:
Is loss/gain of function reflected in RNA-seq transcript counts?
Do LoF/GoF transcripts count toward the RNA-seq TPM count? Or would these LoF/GoF transcripts only be detected by isoform quantification?

Hi Kermit. Can you be more specific about what you mean by "loss-of-function/gain-of-function transcripts"? What scenario are you considering? Are there mutations within the coding regions of the gene, or are there transcript isoforms being expressed that change the function of that RNA? It's tough to know what you're asking about just from what you've written. Feel free to edit the question to add more clarifying information.

@JamesHawley The DNA mutation causes a big enough change that the final protein either no longer works or it takes on a different function

Short answer
Since the effects of any given mutation is so context-dependent (and deciding which effects matter will change between researchers), the process of "variant annotation" is typically separate from RNA-seq mapping/quantification.
If you use a common RNA-seq mapping/quantification tool like STAR, HISAT2, or Kallisto, these counts will combine transcripts with and without the mutation of interest.
This is because for the purpose of mapping, mutations are no different than sequencing errors, and when you're mapping you don't actually know what is what.
Long answer
You might not be able to quantify GoF/LoF at the isoform level, either.
Again, this is because the effects of mutations are so context-dependent.
If you want to know how many transcripts do or don't contain a specific mutation, you'll have to work backwards from that mutation and look at the RNA-seq reads themselves.
You can use tools like Rsamtools, htseq-count, or even IGV for this.
If you want to know whether a given protein's function is affected, you have to consider all the types of mutations than can lead to that protein being affected (in your biological scenario of interest, of course).

Some proteins may lose their binding affinity if a single amino acid is changed, in which case you can look for non-synonymous mutations in the RNA-seq reads.
Some genes may undergo alternative splicing and produce different proteins because of some mutations.
You could find this with a differential isoform expression analysis.
Some genes could be deleted on one allele, leading to a lower expression of that protein.
Biologically, that protein is still expressed, but having less of it makes it less functional in general.
You might be able to find that with copy number estimation or allele-specific expression analysis.

If you're looking at non-coding RNA, what mutations affect their function will be different from what you'd look for if you care about protein function.
All in all, any investigation about what a GoF/LoF mutation looks like and how to go about finding it will depend on you biological scenario.
And you may not find it by just looking at the TPM counts of a gene or its isoforms.


text:
How to chop fasta / bed /peak file genomic segments into smaller fixed or custom genomic intervals
I am trying to split my .bed/.fa file genomic segments into arbitrary smaller overlapping intervals:
Consider a typical line/row in my *.bed file as follows:
chrom. $\ \ \ \ \ $  start $\ \ \ \ \ \ \ \ \ \ \ \ $ end
chr7 $\ \ \ $   127471196 $\ \ \ $    127472363$\ \ \ $      ...
The start and end genomic positions include 1167 base pairs for this line. I need to extend this line, for example, into multiple lines where each line covers exactly 500bp like below:
chr7 $\ \ \ $    127471196 $\ \ \ $    127471696 $\ \ \ $       ...
chr7 $\ \ \ $    127471863  $\ \ \ $   127472363 $\ \ \ $       ...

While each line covers 500pb, between the two new lines they both include the same "range" of genomic positions as the original line including some overlap which is fine (and desired)
Is there an existing tool (perhaps part of bedtools?) that can do something like this?
I know there is the bedops tool that does ALMOST this using the --chop option, where it takes the original segment x1:x2 and chops the original segment into specified fixed bp length intervals, however, if there is any left over towards the end is creates a new non-standard length segment (i.e. x1:x2=10bp and chop length=3, it will give you 4 segments that are length 3, 3, 3 and 1 - what I am looking for is something that adjusts for the last remainder and outputs segments of length 3,3,3 and 3 where the last length 3 segment is right adjusted, while the first three segments are left adjusted or better yet it gives the user the choice of precise redefinition of segment coordinates y1,y2,y3 and y4).  Thank you.

You can run --chop with -x to excise the last piece, adding the --stagger to size the chopped pieces as needed, e.g.:
$ echo -e 'chr1\t0\t10' | bedops --chop 3 --stagger 2 -x -
chr1    0       3
chr1    2       5
chr1    4       7
chr1    6       9

Adjust parameters to taste.
https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#chop-w-chop
https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#stagger-stagger
https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#exclude-x

Thank you @alex - the specific issue with bedops --stagger option is that I have not been able to feed it a negative number, that's what would allow for "overlaps" which is critical!

As shown in the example above, --stagger can give you overlapping elements. However, you also can pipe results to bedops --everything --range N to extend edges by N bases, e.g. ... | bedops --chop X -x - | bedops --everything --range Y - > answer.bed. You'd probably need to use coordinates from fetchChromSizes or similar to trim edge-of-chromosome windows, but that is easy to do via piping to awk or the like.


text:
Relationship between molecular weight and linear peptide
Sorry for this very beginner homework question:
Which of the following linear peptides is consistent with Spectrum = {0 71 99 101 103 128 129 199 200 204 227 230 231 298 303 328 330 332 333}? (Select all that apply.)
CTV

TCE

QCV

AQV

VAQ

TCQ

I tried to use the molecular weight to deduce the answer, but I failed.
In [1]: from rdkit import Chem
   ...: from rdkit.Chem import Descriptors

In [2]: chemical_string = 'A'

In [3]: mol: Chem.Mol = Chem.MolFromFASTA(chemical_string)

In [4]: print( Descriptors.ExactMolWt(mol) )
89.047678464

In [5]: chemical_string = 'C'

In [6]: mol: Chem.Mol = Chem.MolFromFASTA(chemical_string)

In [7]: print( Descriptors.ExactMolWt(mol) )
121.019749464

In [8]: chemical_string = 'E'

In [9]: mol: Chem.Mol = Chem.MolFromFASTA(chemical_string)

In [10]: print( Descriptors.ExactMolWt(mol) )
147.053157768

In [11]: chemical_string = 'T'

In [12]: mol: Chem.Mol = Chem.MolFromFASTA(chemical_string)

In [13]: print( Descriptors.ExactMolWt(mol) )
119.058243148

In [14]: chemical_string = 'Q'

In [15]: mol: Chem.Mol = Chem.MolFromFASTA(chemical_string)

In [16]: print( Descriptors.ExactMolWt(mol) )
146.06914218

In [17]: chemical_string = 'V'

In [18]: mol: Chem.Mol = Chem.MolFromFASTA(chemical_string)

In [19]: print( Descriptors.ExactMolWt(mol) )
117.078978592

In [20]:
```

what kind of spectrum ??

are you trying to accomplish what this server does ? https://iomics.ugent.be/ms2pip/

@pippo1980 Someone told me the following: The first six fragments are the dehydrated mass of each amino acid, so all tripeptides are "consistent". Then it makes no sense again. I think the question might be something like what combination of (cleaved) tripeptides would have dehydrated masses combined to make this set of masses. Then in the dipeptides, the TCE, QCV, VAQ are inconsistent

what does 'the dehydrated mass of each amino acid' means ? AA Mw - H2O Mw ? What about N-term-AA and C-term AA ?

@pippo1980 what do you mean by N-term-AA and C-term AA ?

@pippo1980 See the mass composition analysis

Alanine is 71 = 89 - H2O (18). What about N-term, C-term , first and last aa of polipeptide chain ?

@pippo1980 I am not from chemistry background.  May I ask what N-term, C-term , first and last aa means ?

The free amine end of the chain is called the “N-terminus” or “amino terminus” and the free carboxylic acid end is called the “C-terminus” or “carboxyl terminus, from google could be incorrect

@pippo1980 okay, but how does that help to solve this question ?

Was wondering about your mass composition analysis weight, like N e C term AA weights are missing, unless your machine has a resolution power that cannot distinguish that weights differences


text:
MAGs transcriptomics pipeline question
I have currently a following problem. I have a one sample that I did my metagenomics on (Illumina shotgun + nanopore) and have recovered some high quality MAGs. I also did a metatranscriptomics on the same sample (merging reads with PEAR, removing low quality with Trimmomatic and removing RNA with SortMeRNA). What I would like to do now is to check the expression of the genes in my MAGs with the metatranscriptomics of the sample (hope it's clear enough). I aligned metatranscriptomics reads to each MAG and retained only those reads that aligned (used BOWTIE2). Now I am thinking of doing a functional analysis to check what genes are expressed in each MAG. However, I have no idea how and where to start. Was the alignment even a good idea? How would you tackle this problem?
My idea is to try TransDecoder on these aligned transcriptomics reads to get the ORFs and cluster them with CD-HIT2. Then perform a functional analysis as normal with KEGG/InterProSCAN/Blast etc. Would this be a correct way to go?
Your advice will be highly appreciated.
B


text:
What codes represent what genes?
I want to run experiments on the data used in PatternMarkers & GWCoGAPS for novel data-driven biomarkers via whole transcriptome NMF link
So far, the paper has  reduced the dimension to

ammon's horn
amygdala
anterior cingulate cortex
caudate nucleus
dorsolateral prefrontal cortex
first cervical spinal cord segment
frontal cortex
hypothalamus
putamen
substantia nigra

on the samples/patients GTEX-[11GSP, 11O72, 11OF3, 11TUW, Q2AG, QDT8, XAB4]
I assume they are using it.
Gene read count.    GTEx_Analysis_v6_RNA-seq_RNA-SeQCv1.1.8_gene_reads.gct.gz

A small amount of the data is
Name               Description  GTEX-111CU-1826-SM-5GZYN  GTEX-111FC-0226-SM-5N9B8  GTEX-111VG-2326-SM-5N9BK  ...
ENSG00000223972.4  DDX11L1      0                         0                         0                         ...

From what I can gather, ENSG00000223972.4  DDX11L1  is the DEAD/H-box helicase 11 like 1
Therefore, I assumed that looking for the term first cervical spinal cord segment would give me the name and description.
Therefore, I am wondering if I have the correct data set.

It’s not clear what you’re asking. ENSG00000223972 Is indeed a gene ID. But the terms you list are brain regions, not genes, presumably indicating the sample source


text:
How to merge multiple list into a single dataframe replacing common values with 1 -python
I have multiple lists on python:
a = ["house","garden", "living room", "dog","cat"]

b= ["cat","dog", "chicken"]

c=["house", "garden","bathroom"]

I'd like to create the following dataframe:

a
b
c

house
0
1

garden
0
1

living room
0
0

dog
1
0

cat
1
0

Thanks for the help!

Hi @Marco if you do find the answer below okay (it does work), could you very kindly mark the post as 'accepted'. There's a little grey 'tick' sign to the left of the answer, click that and it goes green. I noted on your past questions you forgotten to do that. Its a small thing but helps everyone, which includes you (you will get +2).

Here  you go, just add column 'C' and you are sorted. Just to note, you must be seeking a pandas solution because merge is a pandas command and you're asking for a Python solution.
import pandas as pd

a = ["house","garden", "living room", "dog","cat"]
b= ["cat","dog", "chicken"]
df = pd.DataFrame(a, columns = ['a'])
df2 = pd.DataFrame(b, columns = ['b'])
dfa = df['a'].value_counts()
dfa.columns = ['a']
dfb = df2['b'].value_counts()
dfb.columns = ['b']
dfa = dfa.to_frame()
dfb = dfb.to_frame()
df3 = dfa.join(dfb).replace(np.nan, 0).astype(int)
print (df3)

Output
             a  b
house        1  0
garden       1  0
living room  1  0
dog          1  1
cat          1  1

and if you want to remove the 'a' column
print (df3.drop('a', axis=1))

             b
house        0
garden       0
living room  0
dog          1
cat          1

Notes The key command here is value_counts and makes a frequency plot. Its output is a Series rather than a DataFrame so it needs converting.
If you want to keep all inputs the commands is
df3 = pd.concat([dfa, dfb]).replace(np.nan, 0).astype(int)

Alternatively,
a = ["house","garden", "living room", "dog","cat"]
b= ["cat","dog", "chicken"]
df = pd.DataFrame(a).value_counts().to_frame('a')
df2 = pd.DataFrame(b).value_counts().to_frame('b')
df3 = df.join(df2).replace(np.nan, 0).astype(int).rename_axis(None)
print (df3.drop('a', axis=1))

I'm aware that this asks specifically for a python solution, but I'd like to add an R answer as well, just in case someone using R has a similar problem:
> a <- c("house","garden", "living room", "dog","cat")
> b <- c("cat","dog", "chicken")
> c <- c("house", "garden","bathroom")

> (result <- data.frame(row.names=a, a=a %in% a, b=a %in% b, c=a %in% c))
               a     b     c
house       TRUE FALSE  TRUE
garden      TRUE FALSE  TRUE
living room TRUE FALSE FALSE
dog         TRUE  TRUE FALSE
cat         TRUE  TRUE FALSE

or for strictly what you asked for, the TRUE/FALSE values can be encouraged into numbers by adding zero:
> (result <- data.frame(a=a, b=(a %in% b) + 0, c=(a %in% c) + 0))
            a b c
1       house 0 1
2      garden 0 1
3 living room 0 0
4         dog 1 0
5         cat 1 0

Both existing answers seem to work to my eye. An alternative approach:
import pandas as pd

a = ["house","garden", "living room", "dog","cat"]
b = ["cat","dog", "chicken"]
c = ["house", "garden","bathroom", "unique"]

First we generate all options available. This way your 'a' list is not restricting the options. We use sets to get unique values.
options = list(set(a + b + c))
wanted = pd.DataFrame({'options' : options} )

Next we loop over all the lists. We use tuples to assign both name and variable.
for name, _ in [('a',a),('b',b),('c',c)]:
    wanted.loc[wanted['options'].isin(_), name ] = 1

Finally we fill the missing values with zeroes.
wanted = wanted.fillna('0')

In case you have large number of lists, you can iterate over ascii_lowercase letters:
my_tuples = ()
from string import ascii_lowercase

for _ in ascii_lowercase[:3]: #<-magic number for number of lists here
     my_tuples = my_tuples + (str(_), eval(_))

Then change for loop as:
for name, _ in my_tuples:
     ....


text:
How to calculate average BLOSUM62 scores?
I can understand the motive behind the BLOSUM62 matrix, this being a pairwise mutation matrix describing aggregate mutations between the 20 amino acids.
However how would you calculate the average BLOSUM62 score within a hydrophobic amino acids (K, D, E) vs hydrophilic amino acids (I, L, V)?
How about between those two groups?

I've not done this precise calculation. I would use point 3, but its complicated to implement.
Just for information BLOSUM62 is a static matrix 20x20 representing relative amino acid mutations. However, I am not exactly sure how the matrix was calculated.  The algorithm looks quite simple, then a BLOSUM62 score can be obtained by comparing the observed matrix to the static matrix, this would provide a score. The paper is here.
Maximum likelihood approaches require a phylogeny.
There are three possibilities I can think of, the first one wouldn't work under a maximum likelihood approach:

Remove all hydrophobic residues from the alignment, check the BLOSUM62 score, remove all hydrophilic check the score etc... again a maximum likelihood approach is not meaningful in this example, but a BLOSUM62 score might be okay;
Recode the amino acids to generate 6-state amino acid data, or recode hydrophobics and in a separate data set hydrophilics and leave the rest uncoded. I've never done it but at a guess I, L, V could be represented by just one amino acid, e.g. L. If the observed BLOSUM62 matrix algorithm was calculated on the observed data the new amino acid e.g. L would have an average score.
Alternatively, manipulate the BLOSUM62 matrix, however its no longer a BLOSUM62 matrix, but it is highly amenable to maximum likelihood analysis. Thus I <-> V and I <-> L is set to zero (or massively negatively weight them) check likelihood, then hydrophilics are set to zero check likelihood. Do the same with the observed versus BLOSUM62 matrix as the control.

Point 3 what a maximum likelihood approach would do is say how close a given matrix is to observed data set for a given tree under a likelihood probability, the smallest is the best. So if you've got loads of matrices its a really cool approach. There is a membrane protein matrix - which is likely heavily skewed towards hydrophilicity. Basically, I'd test all the matrices - but not recode or delete the data set in any way. The matrix that is closer to the data set using likelihood is the best representation of its hydrophilicity vs. hydrophobicity.
Note Just finally mention you can simply calculate the observed amino acid matrix for that data set (via a program) and then average the amino acids. That might be easier and using parsimony its very quick but you'd need a tree. The underlying metric (or metric behind the metric) in all cases is mutations.


text:
Avoiding false anticorrelations between individual genes per cell in single-cell RNAseq
When analysing a single-cell RNAseq dataset, one question I like to ask is: At a single-cell level, are these two genes correlated (expressed together) or anticorrelated (only one or the other is expressed)? However, I am concerned that - given the low read count per cell - it may falsely appear that two genes are anticorrelated when in fact it is just unlikely that both RNAs are detected at the same time.
This problem can partially be dealt with using an imputation technique such as SAVER, to estimate the true relative expression values of undetected genes. Additionally, a statistical test (e.g. permutation) can be performed to ask whether the incidence of cells with both transcripts detected is higher or lower than would be expected, given the total number of trancripts.
My concern, however, is that there may be additional technical aspects of the sequencing process that will exaggerate the apparent anticorrelation between genes. For example, in a cell where the total number of RNA transcripts at any one time is low, a single transcript may be read multiple times (due to multiple primers binding, resulting in multiple cDNAs).
In addition to this, are there any other technical or analytical aspects of ssRNAdeq likely to lead to false anti-correlations?


text:
Merge/Join one data frame (X) to multiple data frames separately
Let us say that we have some data frames named S1,S2,S3... that I want to join independently by column names to dataframe X. This is how they look like:
S1 <- data.frame(
  values = c (1:5),
  numbers = c(11:15))
S1

S2 <- data.frame(
  values = c (1:5),
  numbers = c(11:15))
S2

S3 <- data.frame(
  values = c (1:5),
  numbers = c(16:20))
S3

#Dataframe X
X <- data.frame(
  values = c (26:30),
  numbers = c(501:505))
X

S1
S2
S3
X

Using rbind, it is really easy to merge/join them according to column names. As below:
S1_X<- rbind(S1,X)

S1_X

S2_X<- rbind(S2,X)

S2_X

S3_X<- rbind(S3,X)

S3_X

Now the problem is that I have 500 plus data frames (S1,S2,S3,S4,S5.....S500) that I want to join to dataframe "X" by any function. For example rbind. Is it possible to use some loop/function for that.
What I have best tried is create a list (nested list of dataframes) that I can later split into individual dataframes:
output = lapply(list(S1=S1,S2=S2,S3=S3),rbind,X)

The output is what I expected. So it works The only problem is, how do I specify S1 to S500?
Thank you in advance!

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

Hi, I did specify that I have 500 data frames named S1,S2,S3,S4.....S500 that I want to join to data frame X which is a dataframe too. If you will run the codes above you will find out that I tried to make it simple. Thank you!

@salman1490 I would recommend deleting this question here and reposting on Stack Overflow. It's not really as bioinformatics question and you'll get a much quicker response there.

You can use For Loop here. Here I have created code using your algorithm, first created a list of merged dataframes and then split the list into separate dataframes.
S1 <- data.frame(values = c (1:5), numbers = c(11:15))
S2 <- data.frame(values = c (1:5), numbers = c(11:15))
S3 <- data.frame(values = c (1:5), numbers = c(16:20))
#Dataframe X
X <- data.frame(values = c (26:30), numbers = c(501:505))

#defining a NULL object
S=NULL
#creating a list of merged dataframes
for (i in 1:3){
S[i] = list(rbind(get(paste("S",i,sep='')),X))
names(S)[i] = paste("S",i,"_X",sep='')}

#extracting dataframes
for(i in 1:3){
assign(paste("S",i,"_X",sep = ""),S[[i]])
}