text:
Generating contigs from illumina Paired End Reads (Metagenome)
I'm a wet-lab biologist desperately trying to get to grips with some Bioinformatics.
I have some shotgun Illumina paired-end data from a gut microbiome project I'm working on. I'm considering how to go about generating contigs from reads. I've been recommended to use SPAdes which I'm reading up on at the moment.
Originally I was thinking of using FLASh however I was told it might not be suitable for metagenomic data.
My question simply is, what makes certain programs suitable for merging paired-end reads from a genome but not a metagenome? I find myself a little spoilt for choice (eg. FLASh/SPAdes/BBmap... )  with bioinformatics tools I'm always concerned I'm picking the wrong one.
Should I ideally be running my data through a couple of tools and comparing them?

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

Is the source material metagenomics or culture isolates?

What is the purpose of your project? You have mentioned, "merging paired-end reads", but it's unclear why this would be helpful. If this is for quantifying species or determining proportions, then you'd be better off with a metagenomic mapping tool (e.g. kraken2).

NB: "merging PE reads" (such as e.g. FLASH does) is not the same as "sequence assembly from PE reads" (what MetaSPAdes does). The latter generates contigs, the former generates (some) merged shotgun library inserts. I think you want sequence assembly, i.e. a tool like MetaSPAdes.
Don't be tricked by the loose language of "merging" that is sometimes applied to the sequence assembly problem. It is much more complicated than that.
If you want to learn more about sequence assembly, I suggest you read the sequence assembly wikipedia page.
I would probably recommend MegaHit over metaSPAdes, because it requires much less memory / runtime, and lower frequency of errors.
But generally speaking, both programs are building different flavors of de Bruijn Graphs from metagenomic reads to generate contigs/scaffolds. This is the accepted method by this point for Illumina metagenomic data, because it is the method that is capable of handling the amount of small pieces of data in a tractable fashion. Alternatives such as overlap-layout-consensus graphs or string graphs just don't work as well for the metagenomic problem.
You are probably looking for methods that talk about complicated, computationally complex assembly algorithms and graphical methods for summarizing lots of short reads. You are not looking for simple read processing/alignment tools like BBMap/BBTools or FLASH. To reiterate, simple read merging tools will not give you assembled contigs.

Thank you, this had helped a lot, I'm definitely struggling with some of the subtleties of data processing. My end goal is to do Taxonomic and Functional analysis in MEGAN, but I've been told that it's only designed to use R1 data. I thought I'd attempt to make the most of the sequencing data that I do have to get into a format suitable for MEGAN. Would using  MegaHit/metaSPAdes mean I wouldn't have to put my data through DIAMOND ?  I'm still very new to this and just making sure I have all the steps correct

I am not sure what you mean by only R1 data. I don't think that assembly will make DIAMOND work any better, though it may be helpful for some MEGAN workflows. For examples of both assembled/raw shotgun data for MEGAN+DIAMOND, I'd suggest this ref: https://currentprotocols.onlinelibrary.wiley.com/doi/10.1002/cpz1.59

Thank you, I'm still very much getting to grips with this.  I originally did 150 paired-end sequencing on my samples. I think DIAMOND/MEGAN only accepts a single input file (What I meant by only using Read1 data, I wasn't being clear there.) I thought trying to merge overlapping reads would give me longer DNA sequences for Taxonomic binning.

@dunc4n I'll first note as a tip for the future that having this intention stated in your question would make it a lot easier to understand what you need. I suppose that it's formally possible that read merging would improve things slightly, in practice I am not sure that it makes any difference. I would follow the example of the reference that I linked in my earlier comment, or really any other paper that uses DIAMOND and MEGAN, to understand what to do to preprocess data for those tools. Is there a reason why this won't work for you?

I'll keep that in mind in the future thanks. I don't see any reason why it would work. What had me confused was when I mentioned my plans I was advised that FLASh wouldn't be suitable to generate contigs from my metagenome samples. I had planned to generate contigs from Read1 + Read2  and potentially cat on the output of the unmerged read by FLASh or align them separately and join them post assembly.   I wonder if when I mentioned it originally They thought I was intending to rely on it for contig generation and assembly.   Thank you again for your input it's been really helpful !

@dunc4n  yeah, so I think that the source of the confusion is that merged reads aren't "contigs". Merged reads may represent "contiguous sequences" but the word "contig" has a rather specific meaning in genomics. I'm still not totally clear on the workflow here, but it seems like the issue is resolved at least.


text:
Tutorial for the WGCNA: changes in heatmap colours
I am trying to reproduce the results of the R tutorial of the WGCNA package.
In section I number 5, when generating the heatmap it is quite similar to the one provided by the pdf but the colours of the heatmap differ from the document. Is it possible to change it?
I attach the image together with the code used.
[
Image Obtained

Image of tutorial
nSelect = 400
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select];
sizeGrWindow(9,9)
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")

Could you please provide a link to the WGCNA tutorial?

Please consider the answer by @CorteZero - looks fine to me - and consider upvoting/accepting. Accepting an answer does help with site stats

You can add a color palette to the plot.
e.g.
my_palette <- colorRampPalette(c("red", "orange", "yellow","white"))(n = 100);
TOMplot(plotDiss, selectTree, selectColors, main ="Network heatmap plot, selected genes", col=my_palette)

source: Stack Exchange

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

sort: votes


text:
Downloading genomic protein files from accessions in Python
I am trying to download the _protein.faa.gz files for genomes given their accession numbers through Python. Ideally, I would like to do this without third party libraries. Essentially what I have is a list of just the GCA or GCF accessions. The issue with the ftp site is that it includes the project name along with the accession, which I do not have ahead of time. It would be perfect if I could run downloads from within Python that matched a pattern something like this:
https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/145/295/GCF_000145295.1_*/*_protein.faa.gz
Is it at all possible to do this? Or is it possible to run something similar via efetch in Biopython?
Thanks in advance!

Concept code only
from Bio import Entrez

Entrez.email = "m@M__"
handle = Entrez.efetch(db="protein", id="GCF_000145295", "GCF_000145294", "GCF_000145293", retmode="text",rettype="gb") # check the quotes
records = Entrez.parse(handle)
handle.close()
for record in records:
    print (record)

You can write a loop, flatten it and then insert it to the id= tag.
mylist = ['GCF_000145295', 'GCF_000145294', 'GCF_000145293']
flattened = ','.join(mylist)

The above answer is right - here are some more links that describe what you want to do in more detail:

https://eeob-biodata.github.io/BCB546X-python/06-biopython/
https://medium.com/@insidedna/biopython-to-retrieving-fasta-sequences-from-ncbi-databases-20e35772a25

You're going to want to use the BioPython package as it has functions that do exactly what you're looking for.

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


text:
Extracting clusters, hierarchical clustering - bulk RNASeq
I have been attempting to extract particular clusters, given by hierarchical clustering outputs from the pheatmap function (in R).
Please find the dput output of the dataset in question below:
structure(list(`72v24.logFC` = c(-1.47511280138092, -1.10386310802134,
-1.14476608274903, -1.63563574551239, -0.895377907265352, -1.35216488439027,
-2.04949381135359, -1.0344565728341, -0.967008696515716, -1.11233305998901,
-2.02663681865178, -1.74756562789736, -1.02076838668592, -1.27981009782201,
-1.13181804725663, -1.19537930856639, -2.23476188503608, -1.67500131867632,
-0.950059974858215, -2.3124598902615, -1.93757457440581, -3.24602437386808,
-0.839869104053795, -1.40329424193006, -1.63484798229567, -0.912148120310492,
-1.33467209346367, -1.27031027753064, -0.925498608058671, -0.802956638081314,
-0.85457796397113, -1.0591348681193, -0.825278459132589, -0.671570119524199,
-1.8439768319437, -0.947033608766851, -1.63849733567277, -1.40584263305038,
-1.1317958487115, -1.40724695083764, -1.21842141154588, -2.29248800050516,
-1.41699787526721, -0.879218050664484, -2.25064962332145, 1.4685855647974,
-2.02553848518381, 1.32808748053321, -2.46062484389627, -2.45830217362141,
-0.880536960021976, -0.689060222340137, -0.863262774594, -1.64961312646758,
1.10572346456398, -0.891760767302047, -1.22716047163938, -0.680697741299578,
-1.24551344420046, -0.891918251648935, -1.64722787848732, -0.832845262108387,
-0.937129577660045, -0.891684122691785, 1.96452911491088, -1.59420664438654,
-0.921295293242016, -1.8304480291305, -0.632859040349798, 1.69982819847455,
-0.982668795391973, -1.61816842714933, -1.23602172818196, -2.15356245206668,
-1.67961932563564, 1.15478114880933, 1.32665444151868, -1.04590896017973,
-2.29770026795487, -1.27498150009975, -1.36310217581123, -1.09169550738617,
1.5349413586602, -1.66842984591719, -0.760122936706606, -1.10054393155051,
-1.00344072779106, -0.993259123487476, -1.43136824963825, -1.07214189069438,
-0.990813303020613, -1.51324812740189, -1.23465012540084, -1.00098231964124,
-1.92195764952845, -1.00139179559519, -0.940118366042393, -2.35969198434574,
-0.574626090826994, -0.884870428965888, -0.87601460107501, 1.61336157755029,
-0.777670642457121, -0.7312214522498, -0.812926777785119, -1.13443421373561,
-1.56498420350444, -0.69126972981765, -1.44351026281231, -1.21962275860236,
-2.21155555407497, -1.02015250027788, -0.604620179871691, -1.53013572266942,
-0.790050747397837, -1.17950736026735, -0.781075150161975, -0.888554144538286,
3.59172865862026, -0.8665071699621, -2.97701792242157, -1.44424941630738,
-1.50528882788757, -0.772465957539335, -1.57664185213987, -1.02371443439054,
-0.931807409493221, 1.77728969821973, -1.05455101137715, -1.09765136526264,
-0.881436416213335, -0.985605985458853, -1.31644334991692, 2.26917914828263,
-1.43844640540614, -1.75106843701475, -1.22483991752173, -1.65032906954133,
-0.686319914242903, 2.45473055753775, -1.15594379837188, 1.46397122502826,
-0.91260428248088, -0.726540148411507, -1.18227838019011, -0.940722690371935,
-1.28074043595918, -0.629033607437977, 1.01053893402133, -0.735590704487842,
-2.0620695625823, -0.692251638781327, -0.717818968339447, -0.564023101160995,
1.51541391124355, -1.89444914881727, -0.737931396419173, -1.17763728235316,
-0.887605568116903, -2.8392997898765, -1.48018352596355, 1.60698486226092,
-1.96804216410227, -0.881437183709851, 3.08730764755417, -0.600282670155725,
-0.905719394102218, -1.05553935946213, -1.13508709782814, -0.613381276596874,
-0.642647139688194, -0.861929974132332, -1.30652597213552, -1.73647794867498,
-0.619909931049252, -0.880145154869759, -0.60084263914673, -1.16825573255754,
-0.942919502394457, -1.3709930239358, -0.80595477526712, -0.628115199442632,
-1.11502998682131, -1.01274335471793, -1.16595378001606, -3.36098944183559,
-1.19617345920841, -2.10811718229285, 3.33752122936, -1.67696701892834,
-1.90639199489612, -0.574737209544065, -1.66853541929234, -0.595093777485911,
-0.7964813653225, -0.791102813817634, -0.601987407016077, -1.52559881980829,
1.15480324851026, -0.830167728660421), `72v0.logFC` = c(-5.10937036648522,
-8.41670756910048, -8.18671783299995, -2.79415692883273, -3.62090016917187,
-8.37523504570539, -3.79100494829738, -7.20659380156151, -5.46135634576091,
3.43493215839382, -1.67650075643956, -2.99391647274475, -5.77357941089076,
-5.22675557853878, 3.26436094268851, -5.98693257837598, -5.85864442290848,
-2.39995052889286, -5.19221784614079, -4.72634343497625, -7.46407809296449,
-3.04588366694082, -3.85014181238163, 4.44756662720351, -4.42031795110716,
-2.58478993868624, -5.46121416816988, -2.09364835701952, -5.42397023596482,
-5.26192471384563, -4.48979100280513, -5.14171942073048, -5.5267515390142,
3.88762033289334, -5.88489034628499, -2.85490179670763, -5.31586479048071,
-2.64877712735579, -6.24766069294439, -5.35826818217055, -4.97204642950896,
-2.82989963758335, -4.75016577988789, -3.51684070536132, -3.20931950649022,
-3.62394584217955, -2.81449967053025, -5.68308664428684, -3.11500402887324,
-4.08965188145047, -5.43176189260205, -4.30737780091769, 2.52666133143524,
-3.88213803962076, -5.517506523237, -3.27502948978619, -5.37724717693793,
2.77090845859761, -4.22554364178978, -5.81944792042417, -2.23700806354102,
-1.81804660477942, -3.45613077773548, -4.88444075082709, -5.48877199958173,
-4.90077849970695, -3.79304730634383, -0.953597409660203, -5.573920008042,
-1.83250444732307, -1.73248469438834, -2.35460087734472, -4.43029236446579,
-3.37654864236982, -3.63410171518379, -2.78276966423079, -2.01664452557169,
-3.54834187597395, -2.04395822264753, -2.15693095907293, -4.04759665569609,
-4.87928960250415, 1.1358557257212, -4.89657506999093, -2.42698466771808,
-4.55703295523896, -2.53230009478764, -1.74571568381108, 1.93875826054123,
-2.29178704384166, -1.88631492984573, -4.17811684330152, -4.55155213716622,
1.93093495783903, -3.66429213603731, 2.73817157445051, -3.96977406502995,
-4.18487497756956, 2.17512996773018, -4.41168757499049, -2.96267804428878,
-4.60683474159913, -3.0573109426718, -1.59500391193586, -4.65882130908176,
-2.68903554919404, -4.44676809434802, 1.60275224879078, -3.83760729167497,
-3.09907199540008, -2.37905633686461, -4.29650784478684, -1.88615561094275,
-4.44026540480975, -2.01996858091239, -2.46842498723876, -4.24053886261736,
-4.32561109033911, -2.96097578889787, -3.76708395833968, -2.45229520135946,
-2.33856203921674, -4.23241140779399, -3.48627475173498, 2.25888018977625,
-3.85271202182208, -4.61317704601435, -1.95186634348091, 1.3205477461731,
-2.51561846705872, -1.7031812569795, -4.44717736673796, -3.88478215160905,
-3.27547652595564, -2.81690169904308, -4.83981492357284, -2.4486942125238,
-4.31360736969916, -4.22122930700826, -3.49045260391997, -2.87862947524065,
-4.14268228866617, -2.62252967420329, -4.23902905919912, -4.14777621540367,
-2.86056694414781, -1.70668584273452, 1.3624660569909, -2.59153842674144,
-4.12393693326418, -2.12959469550212, 1.46923996462969, -1.98633421258236,
3.22176343046404, -3.52110060349997, -2.69900299957347, -4.20126962458439,
-4.68116989520068, 1.57592717464553, -3.93999603065316, -2.79446235029201,
1.53128139579193, 2.65531461899762, -1.33692502254506, -1.40309746138702,
7.1503230589606, -3.03232696595877, -3.99799016693196, 1.36106448869886,
2.48321747034381, -1.81760324077147, -2.93792828398095, -4.24318581706171,
-5.02935968459404, -2.8067837982494, -3.18182441332852, -3.17383228514224,
-3.75210155230786, -3.87059449392306, -3.38923922246314, -3.8492188722049,
-4.54207754798088, -3.65594980956956, -4.083692310451, 1.30846845334021,
-2.14067263149232, 0.87594188820449, -2.36549598411797, -1.6026649015159,
3.06122694909985, -2.80125769402892, -1.21061533381979, 2.45621751938726,
-1.29479960958983, -1.32215076834374, -3.40665434275434, -3.87581533996246,
2.25817917663687, -4.47531796954113, -3.0627832955887), `24v0.logFC` = c(-5.63038870666508,
-6.44381761773614, -9.71855520117086, 2.83617729007977, -3.39220401105853,
-9.03206578896647, 4.72316302034493, 2.69600384046815, -3.70005627557764,
1.57004091545159, -2.69837465613229, -1.96398920017231, -6.78881951127235,
2.88395305274553, 2.2815237548239, 4.26299665984366, 1.80759581466942,
-6.13159865127454, 2.31984691801986, 2.42017005000243, 1.90522291464817,
-6.08799743507621, 1.93712248733409, -2.85900752558212, -2.53314737897312,
-5.39825301043452, -2.67723981805717, 1.99371094854688, -6.30018858053764,
-5.51269257077614, 1.22241656195939, 2.21258569930309, -3.34360626488873,
2.84622416535537, 1.62066860166637, -4.78372642162402, -5.70058663575352,
-5.38484766744183, -5.09622306588101, 1.4720608195385, 2.81999896316429,
-5.00082275569591, -4.22688753149231, 3.13575564914666, -3.65108053066376,
-4.96357882349085, 1.41296031244116, 2.02777799934303, -2.69828900301196,
3.3280993785528, -4.91524305899842, 1.33559456145784, -3.20704073441948,
3.4589509748052, -4.76220190139647, -5.41080867153898, -2.52368474702804,
2.35227625126478, -1.66843249760588, -5.98409816921998, 1.97783827997787,
-5.05713668502354, -5.66050145489872, -3.23578154397873, -5.47590590732768,
1.79478248294125, 1.53715323981994, 1.17924158724962, -3.38135082340585,
-6.0533528170201, 3.64133939415158, -5.18911327526881, 1.12836935335145,
-5.50962223997537, -5.13018380364927, 2.43347308480167, 3.09821960202745,
1.08371028206663, 1.15073045150599, 1.66618928654645, -4.57926216582183,
-4.39710926752864, 2.29767598663356, 1.11369483753166, 1.90655785334723,
-4.57402537033236, -1.64795888694238, -4.90903507231585, -4.69225607476014,
-3.78658812613927, 2.71158710397443, 2.0014280958994, 1.42976609883019,
1.69469179645011, 2.02542067939079, 1.29620189704961, -2.11262483065018,
-4.45917113273047, -5.00209626672923, -2.22313144605634, 2.28015100958449,
2.99018492359164, -4.55487718731015, -4.4548605761304, -2.74564764773372,
1.50453692418261, -1.80916071906589, -3.78464149018, 5.44303553967597,
1.63876298449708, -2.55262897048338, -4.52516489294599, 1.86420921266273,
1.92394285952202, -1.87942876800687, -4.72261410600152, -2.57805004195914,
-4.47016439150427, -2.79823793519691, 1.49664827810044, -5.12377787039451,
-1.64302437544027, -2.32977998370269, -1.54481646746192, 1.51575792656483,
-1.66914425228449, -5.67499677983564, -2.41648326328024, 1.4178331810756,
1.34477707145884, 1.6168542644666, 1.51101526992765, 2.85584943327062,
-4.19149110174371, 1.69988651757158, 2.42852927065545, 1.66468253437414,
3.21304236161755, 2.52923647014318, -3.32575501787469, 1.32217968706286,
-4.3404208961354, 1.60065358035513, 1.11559852847827, 2.24622871386089,
-2.07936338453629, -4.096641542765, -1.81072034410725, -1.71317603678557,
3.48120770906848, -3.87801910383696, -4.08020770103266, 1.11478378387699,
-3.46079055745071, 1.7322630648363, 1.54861321913459, 1.96850013686628,
2.50282552968501, -1.81077760463368, 2.02487013804019, 1.95047564703855,
-2.99930567848563, -1.79193606152703, -2.54724389366544, -4.15112678166129,
1.13069494609896, -1.43241741260756, 2.36383913233529, 1.65786717131242,
-2.29507370488367, 1.5786351022915, -3.31042999789919, 2.582031624397,
-3.62438687496853, -4.20056844288066, 1.59787844454997, -4.39077226372811,
-4.50758458577604, -3.5310616440773, 1.34960913912344, 1.56097908933761,
1.80910781056029, -4.1248594098749, 1.5318048796277, -3.97753680837585,
-2.80501778688795, -3.20129716412254, -1.98305552020422, -4.25045728757901,
-4.08742040090535, -4.31190449886363, 1.39895836419766, 1.67102855271528,
-5.13520809951127, -4.31701572736756, 0.832778953035371, 2.30899115735356,
-3.18620970041199, -3.80052092618969, 1.59131320240714)), row.names = c("ENSMUSG00000097417.1",
"ENSMUSG00000084950.7", "ENSMUSG00000021824.12", "ENSMUSG00000085505.1",
"ENSMUSG00000018381.15", "ENSMUSG00000098716.1", "ENSMUSG00000092134.1",
"ENSMUSG00000097272.1", "ENSMUSG00000039206.13", "ENSMUSG00000097860.1",
"ENSMUSG00000047880.14", "ENSMUSG00000087700.1", "ENSMUSG00000094127.1",
"ENSMUSG00000101162.2", "ENSMUSG00000079179.9", "ENSMUSG00000107317.1",
"ENSMUSG00000102451.1", "ENSMUSG00000098609.1", "ENSMUSG00000109669.1",
"ENSMUSG00000086836.1", "ENSMUSG00000092545.7", "ENSMUSG00000106093.1",
"ENSMUSG00000090625.1", "ENSMUSG00000009585.17", "ENSMUSG00000047988.1",
"ENSMUSG00000027950.7", "ENSMUSG00000086430.1", "ENSMUSG00000105941.1",
"ENSMUSG00000004071.12", "ENSMUSG00000104208.1", "ENSMUSG00000022223.8",
"ENSMUSG00000031145.15", "ENSMUSG00000047369.15", "ENSMUSG00000097347.1",
"ENSMUSG00000086680.2", "ENSMUSG00000097660.1", "ENSMUSG00000002109.14",
"ENSMUSG00000108524.1", "ENSMUSG00000067786.16", "ENSMUSG00000021485.13",
"ENSMUSG00000097768.1", "ENSMUSG00000085437.8", "ENSMUSG00000086220.1",
"ENSMUSG00000075576.12", "ENSMUSG00000102299.1", "ENSMUSG00000079293.11",
"ENSMUSG00000090952.2", "ENSMUSG00000040723.14", "ENSMUSG00000027329.9",
"ENSMUSG00000106005.1", "ENSMUSG00000085425.1", "ENSMUSG00000046027.17",
"ENSMUSG00000006517.5", "ENSMUSG00000004221.16", "ENSMUSG00000049775.16",
"ENSMUSG00000109881.1", "ENSMUSG00000023988.8", "ENSMUSG00000092371.1",
"ENSMUSG00000087611.1", "ENSMUSG00000102349.1", "ENSMUSG00000092103.2",
"ENSMUSG00000089636.1", "ENSMUSG00000086360.1", "ENSMUSG00000036138.16",
"ENSMUSG00000045382.6", "ENSMUSG00000090081.1", "ENSMUSG00000107999.1",
"ENSMUSG00000041774.14", "ENSMUSG00000030720.16", "ENSMUSG00000027398.13",
"ENSMUSG00000090246.1", "ENSMUSG00000104823.1", "ENSMUSG00000032528.4",
"ENSMUSG00000100783.1", "ENSMUSG00000082848.1", "ENSMUSG00000021250.13",
"ENSMUSG00000009292.17", "ENSMUSG00000109491.1", "ENSMUSG00000090252.1",
"ENSMUSG00000087443.1", "ENSMUSG00000110148.1", "ENSMUSG00000085754.7",
"ENSMUSG00000064339.1", "ENSMUSG00000091160.1", "ENSMUSG00000108332.1",
"ENSMUSG00000108321.1", "ENSMUSG00000048721.2", "ENSMUSG00000032363.15",
"ENSMUSG00000003868.14", "ENSMUSG00000107075.1", "ENSMUSG00000102961.1",
"ENSMUSG00000096957.1", "ENSMUSG00000093006.2", "ENSMUSG00000109754.1",
"ENSMUSG00000027774.16", "ENSMUSG00000090016.1", "ENSMUSG00000012126.16",
"ENSMUSG00000074136.9", "ENSMUSG00000089875.1", "ENSMUSG00000092190.1",
"ENSMUSG00000074213.2", "ENSMUSG00000086564.9", "ENSMUSG00000055725.11",
"ENSMUSG00000028451.12", "ENSMUSG00000078786.9", "ENSMUSG00000092655.1",
"ENSMUSG00000098482.1", "ENSMUSG00000087442.1", "ENSMUSG00000108574.1",
"ENSMUSG00000087669.1", "ENSMUSG00000053442.4", "ENSMUSG00000101674.6",
"ENSMUSG00000078653.4", "ENSMUSG00000108353.1", "ENSMUSG00000097748.1",
"ENSMUSG00000050010.8", "ENSMUSG00000029623.14", "ENSMUSG00000097060.1",
"ENSMUSG00000050627.13", "ENSMUSG00000092386.1", "ENSMUSG00000085471.1",
"ENSMUSG00000087607.1", "ENSMUSG00000061607.14", "ENSMUSG00000026729.9",
"ENSMUSG00000097854.1", "ENSMUSG00000064330.9", "ENSMUSG00000039753.16",
"ENSMUSG00000027706.8", "ENSMUSG00000029360.3", "ENSMUSG00000103864.1",
"ENSMUSG00000085208.3", "ENSMUSG00000073628.5", "ENSMUSG00000102776.1",
"ENSMUSG00000032487.8", "ENSMUSG00000086998.1", "ENSMUSG00000006728.12",
"ENSMUSG00000021273.9", "ENSMUSG00000062591.5", "ENSMUSG00000040936.14",
"ENSMUSG00000002985.16", "ENSMUSG00000046145.7", "ENSMUSG00000038274.11",
"ENSMUSG00000087366.7", "ENSMUSG00000039452.5", "ENSMUSG00000087660.1",
"ENSMUSG00000087206.1", "ENSMUSG00000079355.3", "ENSMUSG00000086693.1",
"ENSMUSG00000030579.10", "ENSMUSG00000084878.1", "ENSMUSG00000108235.1",
"ENSMUSG00000097077.1", "ENSMUSG00000101438.1", "ENSMUSG00000087367.1",
"ENSMUSG00000064351.1", "ENSMUSG00000086312.1", "ENSMUSG00000097221.1",
"ENSMUSG00000107090.1", "ENSMUSG00000037982.16", "ENSMUSG00000095280.1",
"ENSMUSG00000017837.12", "ENSMUSG00000029422.14", "ENSMUSG00000102564.1",
"ENSMUSG00000026749.11", "ENSMUSG00000067288.12", "ENSMUSG00000025353.4",
"ENSMUSG00000109251.1", "ENSMUSG00000086600.8", "ENSMUSG00000052214.9",
"ENSMUSG00000090330.2", "ENSMUSG00000090009.1", "ENSMUSG00000059791.14",
"ENSMUSG00000101695.1", "ENSMUSG00000103698.1", "ENSMUSG00000097472.1",
"ENSMUSG00000084128.10", "ENSMUSG00000026156.8", "ENSMUSG00000075297.10",
"ENSMUSG00000055692.18", "ENSMUSG00000038502.16", "ENSMUSG00000015942.9",
"ENSMUSG00000106478.1", "ENSMUSG00000039307.16", "ENSMUSG00000054556.6",
"ENSMUSG00000045576.16", "ENSMUSG00000037936.15", "ENSMUSG00000087486.1",
"ENSMUSG00000082180.1", "ENSMUSG00000014470.4", "ENSMUSG00000056078.5",
"ENSMUSG00000049928.15", "ENSMUSG00000013997.10", "ENSMUSG00000105198.1",
"ENSMUSG00000086714.1", "ENSMUSG00000109274.1", "ENSMUSG00000100548.1",
"ENSMUSG00000053706.2", "ENSMUSG00000109854.1", "ENSMUSG00000060802.8",
"ENSMUSG00000093598.1"), class = "data.frame")

Please note, I have checked out many posts regarding this, most notably: https://www.biostars.org/p/287512/#287518, which highlights how to perform this procedure. Regardless, there are some unanswered questions.
The problem is the following: After running the code suggested in the post, I do indeed obtain x number of clusters, depending on the height I decide to cut the tree at. However, they are in an incomprehensible order, making it extremely difficult to extract the clusters as they are observed in the heatmap. Please note, the dataframe I have uploaded is only a snippet of the true data - and is called logFC_mat, in the code below:
pheatmap(logFC_mat,
         show_rownames=F, cluster_cols=F, cluster_rows=T, scale="row",
         cex=1, clustering_distance_rows="euclidean", cex=1,
         clustering_distance_cols="euclidean",
         clustering_method="complete", border_color=FALSE) -> heatof

#make clusters based on h=1.5
sort(cutree(heatof$tree_row, h=1.5)) -> hi_clusters

#check if clusters are in same order as in heatmap
hi_clusters[hi_clusters==1] %>% names() -> clust1
hi_clusters[hi_clusters==2] %>% names() -> clust2

length(clust1) #56
length(clust2) #78
all(clust1 %in%
      rownames(logFC_mat[heatof$tree_row$order[1:56],])) #TRUE

all(clust2 %in%
      rownames(logFC_mat[heatof$tree_row$order[57:134],])) #FALSE

Now, from visual inspection of the heatmap, it becomes clear that cluster 2 is in fact the last 78 genes rather than the middle cluster. The problem is, why is it labeled as cluster 2 then !? In addition, this is a snippet of the data. The real data contains in fact thousands of rows, and so such visual deductions can simply be impossible.
Now for the question: How to extract the clusters, such that they are in the SAME ORDER, as what is seen in the heatmap (Top to bottom)?
Thanks.

I would recommend not to use h to cut the tree but to use k, as k instructs the number of clusters to obtain, and then leaves it to the function to find the correct height for that. From there on you can easily extract the genes per cluster, for example with this code obtaining the first three clusters:

heatof <-
pheatmap(logFC_mat,
         show_rownames=F, cluster_cols=F, cluster_rows=T, scale="row",
         cex=1, clustering_distance_rows="euclidean", cex=1,
         clustering_distance_cols="euclidean",
         clustering_method="complete", border_color=FALSE)

#/ extract the top three clusters
k <- 3
sort(cutree(heatof$tree_row, k=3)) -> hi_clusters

genes_per_cluster <- do.call(rbind, lapply(unique(hi_clusters), function(x){

  data.frame(cluster=x, gene=names(hi_clusters[hi_clusters %in% x]))

}))

genes_per_cluster
#>     cluster                  gene
#> 1         1  ENSMUSG00000097417.1
#> 2         1  ENSMUSG00000084950.7
#> 3         1 ENSMUSG00000021824.12
#> 4         1 ENSMUSG00000018381.15
#> 5         1  ENSMUSG00000098716.1
#> 6         1 ENSMUSG00000039206.13
#> 7         1 ENSMUSG00000047880.14
#> 8         1  ENSMUSG00000094127.1
#> 9         1  ENSMUSG00000098609.1
#> 10        1  ENSMUSG00000106093.1
#> 11        1 ENSMUSG00000009585.17
#> 12        1  ENSMUSG00000027950.7
#> 13        1 ENSMUSG00000004071.12
#> 14        1  ENSMUSG00000104208.1
#> 15        1 ENSMUSG00000047369.15
#> 16        1  ENSMUSG00000097660.1
#> 17        1 ENSMUSG00000002109.14
#> 18        1  ENSMUSG00000108524.1
#> 19        1 ENSMUSG00000067786.16
#> 20        1  ENSMUSG00000085437.8
#> 21        1  ENSMUSG00000086220.1
#> 22        1  ENSMUSG00000102299.1
#> 23        1 ENSMUSG00000079293.11
#> 24        1  ENSMUSG00000085425.1
#> 25        1  ENSMUSG00000006517.5
#> 26        1 ENSMUSG00000049775.16
#> 27        1  ENSMUSG00000109881.1
#> 28        1  ENSMUSG00000102349.1
#> 29        1  ENSMUSG00000089636.1
#> 30        1  ENSMUSG00000086360.1
#> 31        1 ENSMUSG00000036138.16
#> 32        1  ENSMUSG00000045382.6
#> 33        1 ENSMUSG00000030720.16
#> 34        1 ENSMUSG00000027398.13
#> 35        1  ENSMUSG00000104823.1
#> 36        1  ENSMUSG00000100783.1
#> 37        1  ENSMUSG00000082848.1
#> 38        1  ENSMUSG00000110148.1
#> 39        1  ENSMUSG00000085754.7
#> 40        1  ENSMUSG00000108321.1
#> 41        1 ENSMUSG00000032363.15
#> 42        1 ENSMUSG00000003868.14
#> 43        1  ENSMUSG00000107075.1
#> 44        1  ENSMUSG00000074136.9
#> 45        1  ENSMUSG00000089875.1
#> 46        1 ENSMUSG00000055725.11
#> 47        1 ENSMUSG00000028451.12
#> 48        1  ENSMUSG00000078786.9
#> 49        1  ENSMUSG00000087442.1
#> 50        1  ENSMUSG00000053442.4
#> 51        1  ENSMUSG00000101674.6
#> 52        1  ENSMUSG00000097748.1
#> 53        1  ENSMUSG00000050010.8
#> 54        1 ENSMUSG00000029623.14
#> 55        1  ENSMUSG00000097060.1
#> 56        1 ENSMUSG00000050627.13
#> 57        1  ENSMUSG00000085471.1
#> 58        1 ENSMUSG00000039753.16
#> 59        1  ENSMUSG00000027706.8
#> 60        1  ENSMUSG00000032487.8
#> 61        1 ENSMUSG00000002985.16
#> 62        1 ENSMUSG00000038274.11
#> 63        1  ENSMUSG00000087206.1
#> 64        1  ENSMUSG00000079355.3
#> 65        1  ENSMUSG00000086693.1
#> 66        1 ENSMUSG00000030579.10
#> 67        1  ENSMUSG00000108235.1
#> 68        1  ENSMUSG00000097077.1
#> 69        1  ENSMUSG00000087367.1
#> 70        1 ENSMUSG00000037982.16
#> 71        1 ENSMUSG00000029422.14
#> 72        1 ENSMUSG00000026749.11
#> 73        1 ENSMUSG00000067288.12
#> 74        1  ENSMUSG00000090330.2
#> 75        1 ENSMUSG00000059791.14
#> 76        1  ENSMUSG00000103698.1
#> 77        1  ENSMUSG00000097472.1
#> 78        1  ENSMUSG00000026156.8
#> 79        1 ENSMUSG00000075297.10
#> 80        1 ENSMUSG00000055692.18
#> 81        1 ENSMUSG00000039307.16
#> 82        1 ENSMUSG00000045576.16
#> 83        1  ENSMUSG00000087486.1
#> 84        1  ENSMUSG00000014470.4
#> 85        1  ENSMUSG00000056078.5
#> 86        1 ENSMUSG00000049928.15
#> 87        1  ENSMUSG00000086714.1
#> 88        1  ENSMUSG00000109274.1
#> 89        1  ENSMUSG00000109854.1
#> 90        1  ENSMUSG00000060802.8
#> 91        2  ENSMUSG00000085505.1
#> 92        2  ENSMUSG00000092134.1
#> 93        2  ENSMUSG00000097272.1
#> 94        2  ENSMUSG00000087700.1
#> 95        2  ENSMUSG00000101162.2
#> 96        2  ENSMUSG00000107317.1
#> 97        2  ENSMUSG00000102451.1
#> 98        2  ENSMUSG00000109669.1
#> 99        2  ENSMUSG00000086836.1
#> 100       2  ENSMUSG00000092545.7
#> 101       2  ENSMUSG00000090625.1
#> 102       2  ENSMUSG00000047988.1
#> 103       2  ENSMUSG00000086430.1
#> 104       2  ENSMUSG00000105941.1
#> 105       2  ENSMUSG00000022223.8
#> 106       2 ENSMUSG00000031145.15
#> 107       2  ENSMUSG00000086680.2
#> 108       2 ENSMUSG00000021485.13
#> 109       2  ENSMUSG00000097768.1
#> 110       2 ENSMUSG00000075576.12
#> 111       2  ENSMUSG00000090952.2
#> 112       2 ENSMUSG00000040723.14
#> 113       2  ENSMUSG00000027329.9
#> 114       2  ENSMUSG00000106005.1
#> 115       2 ENSMUSG00000046027.17
#> 116       2 ENSMUSG00000004221.16
#> 117       2  ENSMUSG00000023988.8
#> 118       2  ENSMUSG00000087611.1
#> 119       2  ENSMUSG00000092103.2
#> 120       2  ENSMUSG00000090081.1
#> 121       2  ENSMUSG00000107999.1
#> 122       2 ENSMUSG00000041774.14
#> 123       2  ENSMUSG00000090246.1
#> 124       2  ENSMUSG00000032528.4
#> 125       2 ENSMUSG00000021250.13
#> 126       2 ENSMUSG00000009292.17
#> 127       2  ENSMUSG00000109491.1
#> 128       2  ENSMUSG00000090252.1
#> 129       2  ENSMUSG00000087443.1
#> 130       2  ENSMUSG00000064339.1
#> 131       2  ENSMUSG00000091160.1
#> 132       2  ENSMUSG00000108332.1
#> 133       2  ENSMUSG00000048721.2
#> 134       2  ENSMUSG00000102961.1
#> 135       2  ENSMUSG00000096957.1
#> 136       2  ENSMUSG00000093006.2
#> 137       2 ENSMUSG00000027774.16
#> 138       2 ENSMUSG00000012126.16
#> 139       2  ENSMUSG00000092190.1
#> 140       2  ENSMUSG00000074213.2
#> 141       2  ENSMUSG00000086564.9
#> 142       2  ENSMUSG00000092655.1
#> 143       2  ENSMUSG00000098482.1
#> 144       2  ENSMUSG00000108574.1
#> 145       2  ENSMUSG00000087669.1
#> 146       2  ENSMUSG00000078653.4
#> 147       2  ENSMUSG00000108353.1
#> 148       2  ENSMUSG00000092386.1
#> 149       2  ENSMUSG00000087607.1
#> 150       2 ENSMUSG00000061607.14
#> 151       2  ENSMUSG00000026729.9
#> 152       2  ENSMUSG00000064330.9
#> 153       2  ENSMUSG00000103864.1
#> 154       2  ENSMUSG00000085208.3
#> 155       2  ENSMUSG00000073628.5
#> 156       2  ENSMUSG00000102776.1
#> 157       2  ENSMUSG00000086998.1
#> 158       2 ENSMUSG00000006728.12
#> 159       2  ENSMUSG00000021273.9
#> 160       2  ENSMUSG00000062591.5
#> 161       2 ENSMUSG00000040936.14
#> 162       2  ENSMUSG00000046145.7
#> 163       2  ENSMUSG00000087366.7
#> 164       2  ENSMUSG00000039452.5
#> 165       2  ENSMUSG00000087660.1
#> 166       2  ENSMUSG00000084878.1
#> 167       2  ENSMUSG00000101438.1
#> 168       2  ENSMUSG00000064351.1
#> 169       2  ENSMUSG00000086312.1
#> 170       2  ENSMUSG00000097221.1
#> 171       2  ENSMUSG00000107090.1
#> 172       2  ENSMUSG00000095280.1
#> 173       2 ENSMUSG00000017837.12
#> 174       2  ENSMUSG00000109251.1
#> 175       2  ENSMUSG00000086600.8
#> 176       2  ENSMUSG00000090009.1
#> 177       2  ENSMUSG00000101695.1
#> 178       2 ENSMUSG00000084128.10
#> 179       2 ENSMUSG00000038502.16
#> 180       2  ENSMUSG00000015942.9
#> 181       2  ENSMUSG00000106478.1
#> 182       2  ENSMUSG00000054556.6
#> 183       2  ENSMUSG00000082180.1
#> 184       2 ENSMUSG00000013997.10
#> 185       2  ENSMUSG00000100548.1
#> 186       2  ENSMUSG00000053706.2
#> 187       2  ENSMUSG00000093598.1
#> 188       3  ENSMUSG00000097860.1
#> 189       3  ENSMUSG00000079179.9
#> 190       3  ENSMUSG00000097347.1
#> 191       3  ENSMUSG00000092371.1
#> 192       3  ENSMUSG00000109754.1
#> 193       3  ENSMUSG00000090016.1
#> 194       3  ENSMUSG00000097854.1
#> 195       3  ENSMUSG00000029360.3
#> 196       3  ENSMUSG00000102564.1
#> 197       3  ENSMUSG00000025353.4
#> 198       3  ENSMUSG00000052214.9
#> 199       3 ENSMUSG00000037936.15
#> 200       3  ENSMUSG00000105198.1

Created on 2022-07-07 by the reprex package (v2.0.1)

```

Thanks for the answer. I have no problem retrieving the genes in the clusters, just identifying the genes associated with the clusters observed in the heatmap. E.g. in my example, the first cluster (cluster 1) is the first cluster in the heatmap (on the top), however my second cluster is not my second cluster (cluster 2)...


text:
Different line length in fasta file
I am currently using VEP for variant annotation.
I am facing an error as below:
[E::fai_build_core] Different line length in sequence 'Pn9'

I understand there is an issue with the difference in line length of Pn9 in the  fasta file. However, the sequences of fasta file all have different line length. I don't get why there is an error specifically on Pn9.
I have tried using both snpEff and Annovar but doesn't work. Any thoughts on this would be really appreciated.
I am attaching the length of sequence here :

Thanks 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.

This could be a bug in the fasta file reader component of fai_build_core. The fasta file format is very relaxed, and allows some odd constructions, for example:
>cb3c6a5f-17dd-4cd6-a9ca-02124d2a32de
AAGAAATTGACCTCCGTAACCCCAAAGGAGACAAAAAGAATGCTGGCTTTACTCCC
CAGGAGCGCCAGGTTTTGGGGAACTGCCTACAAGCTGTACCATTGGCTGACAGCTT
TTAGGGCATCTACCCCCAACACTGCTTACACTTACACTTTCTGGACTTACATG
ATGAATGCCCGCTCTAAGAATGTTGGTTGGCGCCTTGATTACTTTTTGCTTTCC
CACTCTCTTTTTACCTGCATTGTGACAGCAAGATCCGGTCAAGGCTCTTG
GCAGTACACTGTCCCATCACCGCTACCTGTGATACCCCTCCCTGCAGTAGCTT
CCTGCCTGGGAGATGGCTCTCTCTGCAGAAGTCTGGTGTTTTAGCC
TTCAGGTGTTTGGTTTGTATGTGCTCCCTCATTTAAACATTAAACCAAACTT
CTGGTTTCCTTTAGACAATCCAAGAGAAATAAAAGGCCCTAC
TTTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA

A few programs will complain about fasta files that look like this, because it usually means that data has been corrupted. These programs only accept fasta files that have sequences with identical line lengths for all but the last line:
>cb3c6a5f-17dd-4cd6-a9ca-02124d2a32de
AAGAAATTGACCTCCGTAACCCCAAAGGAGACAAAAAGAATGCTGGCTTTACTCCC
CAGGAGCGCCAGGTTTTGGGGAACTGCCTACAAGCTGTACCATTGGCTGACAGCTT
TTAGGGCATCTACCCCCAACACTGCTTACACTTACACTTTCTGGACTTACATGATG
AATGCCCGCTCTAAGAATGTTGGTTGGCGCCTTGATTACTTTTTGCTTTCCCACTC
TCTTTTTACCTGCATTGTGACAGCAAGATCCGGTCAAGGCTCTTGGCAGTACACTG
TCCCATCACCGCTACCTGTGATACCCCTCCCTGCAGTAGCTTCCTGCCTGGGAGAT
GGCTCTCTCTGCAGAAGTCTGGTGTTTTAGCCTTCAGGTGTTTGGTTTGTATGTGC
TCCCTCATTTAAACATTAAACCAAACTTCTGGTTTCCTTTAGACAATCCAAGAGAA
ATAAAAGGCCCTACTTTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA

One way to work around this input error problem is to get rid of all line breaks in the sequence of a fasta file, but then you run the risk of programs not properly being able to read all the characters in a line. A better approach is to feed the fasta file through a less restrictive data cleaning program and convert the lines back to all have a consistent length.
There's a good chance this is not the problem you are experiencing, but without more information about your input data, all we have to go on are the error messages and the sequence lengths [not line lengths] you have provided.

Thank you so much for your reply. I was able to rectify the error. I was doing some unwanted steps while zipping the fasta file.


text:
Get number of reads with a single, (almost) exact match to the full length of a reference sequence
I can't find an answer to this in previous questions, so hoping someone can help me now.
We have nanopore sequenced a PCR product, and I've filtered our reads to +/- 10 bp of the expected product length (~4,400bp), leaving ~55,000 reads.
We now want to work out how many of the reads correspond to each of 6 different potential arrangements. The section we're interested in is in the middle of the 4,400 bp PCR product, and each different potential arrangement is ~1,500 bp. We have reference sequences for each of the 1,500 bp arrangements. Basically, we want to know what portion of our starting cells had each of the different arrangements of genes in this operon (e.g arrangement 1=30%, arrangement 2=15%, etc.).
This seems like it should be a really easy thing to do using read mapping to count alignments for each reference separately, but I just can't get my head around it. No matter what options I've tried with minimap2 and samtools, every read maps every time on every different arrangement. I think the problem is that the genes in our reference sequences are all present in every read, just in different arrangements. So we need to only count alignments which align to the whole 1,500 bp reference sequence, not partial hits or hits with huge gaps in them.
I have tried the minimap2 options -r, -G and --secondary=no to restrict the output to single alignments with no large gaps, and combined that with samtools view -m 1000 to output only alignments longer than 1,000 bp (this is far stricter than it needs to be to only output alignments to any of the different potential arrangements). But still, every time, every read is mapping to every reference, usually multiple times.
What am I doing wrong?! Am I totally over-complicating this, or doing something obviously stupid? Almost feel like using blastn might be better at this point...
TL;DR: how do I use minimap2 and/or samtools (or any other tool better suited to this) to count how many of our ~4,400 bp reads contain only full-length, almost exact matches to our 1,500 bp reference sequence, excluding any partial hits?
My latest attempt (broken into individual steps, though I'd usually pipe the samtools steps together, also not bothered with samtools sort at this point):
minimap2 -ax map-ont -r10 -N5 --secondary=no --no-long-join reference.fasta reads.fastq > reads.minimap.sam

samtools view -m 1000 -b -o reads.minimap.bam reads.minimap.sam

samtools -c reads.minimap.bam

Because minimap2 is a minimiser-based mapper, it's less useful for determining mapping accuracy to single-base precision.
Where precision matters, you can use LAST, bearing in mind that it calculates local alignments, and can report multiple hits per query/reference pair. The usual approach I use for LAST mapping is as follows:

lastdb -uNEAR -R01 reference.fa reference.fa
last-train -Q 1 reference.fa <(cat reads.fastq | head -n 40000) > train.mat
lastal -p train.mat reference.fa reads.fastq > reads_vs_ref.maf

That final MAF file can be inspected for mapping. The default output produced from LAST is not very human-readable, so I've written a perl script that converts MAF results into a CSV format with reference-relative statistics that I have found useful:
maf2csv.pl reads_vs_ref.maf > reads_vs_ref.csv
This can be inspected in a spreadsheet viewer (e.g. LibreOffice Calc, or Excel) to look at the mapped locations, mapped fraction, identities, etc..
Alternatively, the MAF results can be converted to SAM/BAM format using maf-convert (bearing in mind the previously-mentioned caveats about multiple hits):
maf-convert sam reads_vs_ref.maf | samtools view -b -T reference.fa | samtools sort > reads_vs_ref.bam


text:
DNA genome string reconstruction from k-mer
I have the following quiz question, but the Pattern1 for both (ACC|ATA) and (CGA|ACT) are unique (just grep for ACC and CGA), which means it might be difficult to find the starting pattern for the DNA genome string in the following case. Do I need something like GapReduce algorithm ?
Quiz : There is a single (linear) string with the following (3,1)-mer composition. Find it.
(ACC|ATA) (ACT|ATT) (ATA|TGA) (ATT|TGA) (CAC|GAT) (CCG|TAC) (CGA|ACT)

(CTG|AGC) (CTG|TTC) (GAA|CTT) (GAT|CTG) (GAT|CTG) (TAC|GAT) (TCT|AAG)

(TGA|GCT) (TGA|TCT) (TTC|GAA)

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.

As shown in the picture in the question, just align the (k, d)-mers vertically without the help of assembler software, then compare whether they conflict with each other vertically. Maybe the software to do such comparison uses https://en.wikipedia.org/wiki/Smith%E2%80%93Waterman_algorithm


text:
Jitter violin coloured by sample id
I want to have a jitter violin coloured by sample ids but I get this weird plot
> str(long)
Classes ‘data.table’ and 'data.frame':  156 obs. of  3 variables:
 $ gene    : chr  "gene2" "gene3" "gene4" "gene1" ...
 $ variable: chr  "1016" "1016" "1016" "1016" ...
 $ value   : num  5.82 8.41 11.22 7.21 5.61 ...
 - attr(*, ".internal.selfref")=<externalptr>
>

> p <- ggplot(long, aes(x=gene, y=value,colour=variable))
> p + geom_violin(width = 1.3) + geom_quasirandom(alpha = 0.2, width = 0.2)+ theme_classic()
There were 21 warnings (use warnings() to see them)

I can not figure out where I am doing something wrong that I get this instead of getting a proper violin plot by dots coloured with sample IDs

What are the warnings()?

I see 17 different variables and 4 genes on the plot this are 153 values for 1 measurement for each sample and gene but your data.frame has only 156 rows. So at most you have 3 repeated measures for the same gene and variable. You can't make a (meaningful) boxplot of 1 value, so it is expected that very few boxplots can be drawn

What happens if you group by gene?
> long |>
  ggplot() +
  aes(x=gene, y=value, colour=variable, group=gene) +
  geom_violin(width = 1.3) +
  geom_quasirandom(alpha = 0.2, width = 0.2) +
  theme_classic()

@llrs raised an insightful point in their comment. There's a good chance that the "weird plot" that you see is a side effect of trying to display individual "grouped-by-variable" boxplots that have too little data. Without knowing what you would like the boxplot to look like, it's hard to create something to improve the presentation.


text:
AttributeError [in fastStructure's chooseK.py]: module 'vars' has no attribute 'insum'
I've been working with the fastStructure program and am on the step of analyzing model complexity using the provided chooseK.py script. I have been running into the following error:
Traceback (most recent call last):
  File "chooseK.py", line 109, in <module>
    bestKs = parse_varQs(files)
  File "chooseK.py", line 54, in parse_varQs
    Q = Q/utils.insum(Q,[1])
AttributeError: module 'vars' has no attribute 'insum'

I have run into 'AttributeError' messages in the past and have attempted the solution of renaming files in my working directory that may be shadowing the real module (similar to: https://stackoverflow.com/questions/36530726/using-numpy-module-object-has-no-attribute-array), but this hasn't fixed anything. Reinstalling Numpy/Scipy/fastStructure has also been a dead end.
Any suggestions are much appreciated! Happy to provide more info as needed.
-A

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.

did you build the library extensions? The module vars.utils are imported as utils, so the attribute insum is undetected. https://stackoverflow.com/questions/46312470/difference-between-methods-and-attributes-in-python

Could you try loading fastStructure via conda?
conda create --name faststruc python=3.9 # there can be issues with 3.10
conda activate faststruc
conda install -c bioconda faststructure

It will then ask you whether you want to install numpy and scipy automatically. Furthermore all dependencies are installed and connectivity will have already been tested. You can easily delete everything afterwards by deleting the environment faststruc so keeps everything clean.

The code in question is here:
https://programtalk.com/python-examples/vars.utils.insum/
Its a numpy calculation, what they are doing is navigating Python by leveraging the power of numpy, which is very fast.
Code ...
bestKs = []
for file in files:
    handle = open(file,'r')
    Q = np.array([map(float,line.strip().split()) for line in handle])
    Q = Q/utils.insum(Q,[1])
    handle.close()
    N = Q.shape[0]
    C = np.cuemsum(np.sort(Q.sum(0))[::-1])
    bestKs.append(np.sum(C<N-1)+1)

return bestKs

The code looks good, they know their stuff.


text:
R ggtree apply colors to branches with multiple colors
This question was also asked on StackOverflow
I have a phylogenetic tree made with a Newick format
((a:1,b:1):2, (c:1, d:1):3):1;

The output will be

I have drawn this with http://etetoolkit.org/treeview/, but output with R will be same.
What I want to do is apply colors to a certain branch.
That would be

For example,

Since the Newick input doesn't include any information other than the branch (c,d) has length 3,
is there any way I can include multiple color information on the branch?
I usually use ggtree for making trees, but other suggestions can be helpful.
The main script for this kind of tree building is like this, but I can't figure out how this works.
(Full script Github link: https://github.com/seongyeol-park/Human_Lineage_Tracing/blob/main/Rscripts/HSC_commitment.R)
ggtree(tree) %<+% m_lng_dt + theme_tree2()+
  coord_cartesian(xlim = c(-1,35))+
  geom_text2(aes(subset = is.na(Source_class2) == F, x=x-n_pointmt+3, y = y+0.25, label = Source_class2), size=2.5, color="red")+
  geom_segment2(aes(subset = is.na(VAF1) == F, x=x-n_pointmt, xend = x -n_pointmt +1, y = y, yend = y, color = VAF1))+
  geom_text2(aes(subset = is.na(VAF1) == F, x=x-n_pointmt, y = y, label = VAF1), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF2) == F, x=x-n_pointmt+1, xend = x -n_pointmt +2, y = y, yend = y, color = VAF2))+
  geom_text2(aes(subset = is.na(VAF2) == F, x=x-n_pointmt+1, y = y, label = VAF2), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF3) == F, x=x-n_pointmt+2, xend = x -n_pointmt +3, y = y, yend = y, color = VAF3))+
  geom_text2(aes(subset = is.na(VAF3) == F, x=x-n_pointmt+2, y = y, label = VAF3), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF4) == F, x=x-n_pointmt+3, xend = x -n_pointmt +4, y = y, yend = y, color = VAF4))+
  geom_text2(aes(subset = is.na(VAF4) == F, x=x-n_pointmt+3, y = y, label = VAF4), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF5) == F, x=x-n_pointmt+4, xend = x -n_pointmt +5, y = y, yend = y, color = VAF5))+
  geom_text2(aes(subset = is.na(VAF5) == F, x=x-n_pointmt+4, y = y, label = VAF5), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF6) == F, x=x-n_pointmt+5, xend = x -n_pointmt +6, y = y, yend = y, color = VAF6))+
  geom_text2(aes(subset = is.na(VAF6) == F, x=x-n_pointmt+5, y = y, label = VAF6), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF7) == F, x=x-n_pointmt+6, xend = x -n_pointmt +7, y = y, yend = y, color = VAF7))+
  geom_text2(aes(subset = is.na(VAF7) == F, x=x-n_pointmt+6, y = y, label = VAF7), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF8) == F, x=x-n_pointmt+7, xend = x -n_pointmt +8, y = y, yend = y, color = VAF8))+
  geom_text2(aes(subset = is.na(VAF8) == F, x=x-n_pointmt+7, y = y, label = VAF8), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF9) == F, x=x-n_pointmt+8, xend = x -n_pointmt +9, y = y, yend = y, color = VAF9))+
  geom_text2(aes(subset = is.na(VAF9) == F, x=x-n_pointmt+8, y = y, label = VAF9), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF10) == F, x=x-n_pointmt+9, xend = x -n_pointmt +10, y = y, yend = y, color = VAF10))+
  geom_text2(aes(subset = is.na(VAF10) == F, x=x-n_pointmt+9, y = y, label = VAF10), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF11) == F, x=x-n_pointmt+10, xend = x -n_pointmt +11, y = y, yend = y, color = VAF11))+
  geom_text2(aes(subset = is.na(VAF11) == F, x=x-n_pointmt+10, y = y, label = VAF11), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF12) == F, x=x-n_pointmt+11, xend = x -n_pointmt +12, y = y, yend = y, color = VAF12))+
  geom_text2(aes(subset = is.na(VAF12) == F, x=x-n_pointmt+11, y = y, label = VAF12), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF13) == F, x=x-n_pointmt+12, xend = x -n_pointmt +13, y = y, yend = y, color = VAF13))+
  geom_text2(aes(subset = is.na(VAF13) == F, x=x-n_pointmt+12, y = y, label = VAF13), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF14) == F, x=x-n_pointmt+13, xend = x -n_pointmt +14, y = y, yend = y, color = VAF14))+
  geom_text2(aes(subset = is.na(VAF14) == F, x=x-n_pointmt+13, y = y, label = VAF14), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF15) == F, x=x-n_pointmt+14, xend = x -n_pointmt +15, y = y, yend = y, color = VAF15))+
  geom_text2(aes(subset = is.na(VAF15) == F, x=x-n_pointmt+14, y = y, label = VAF15), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF16) == F, x=x-n_pointmt+15, xend = x -n_pointmt +16, y = y, yend = y, color = VAF16))+
  geom_text2(aes(subset = is.na(VAF16) == F, x=x-n_pointmt+15, y = y, label = VAF16), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF17) == F, x=x-n_pointmt+16, xend = x -n_pointmt +17, y = y, yend = y, color = VAF17))+
  geom_text2(aes(subset = is.na(VAF17) == F, x=x-n_pointmt+16, y = y, label = VAF17), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF18) == F, x=x-n_pointmt+17, xend = x -n_pointmt +18, y = y, yend = y, color = VAF18))+
  geom_text2(aes(subset = is.na(VAF18) == F, x=x-n_pointmt+17, y = y, label = VAF18), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF19) == F, x=x-n_pointmt+18, xend = x -n_pointmt +19, y = y, yend = y, color = VAF19))+
  geom_text2(aes(subset = is.na(VAF19) == F, x=x-n_pointmt+18, y = y, label = VAF19), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF20) == F, x=x-n_pointmt+19, xend = x -n_pointmt +20, y = y, yend = y, color = VAF20))+
  geom_text2(aes(subset = is.na(VAF20) == F, x=x-n_pointmt+19, y = y, label = VAF20), size=2.5)+
  scale_color_gradientn(colors = my_palette, breaks = col_breaks)+
  theme(axis.text = element_text(size=15, angle=45, hjust=1), plot.title = element_text(size=20), plot.margin= unit(c(0.2,0,0.2,0.2),"cm"))

I added a random colour tree, which is what I think you might want excepting the legend.

The simple solution is Figtree a GUI that easily colours trees on a 'point and click' basis. In linux this can be driven via CLI, which is useful.
The easy coding solution is to use extended Newick format (NHX) to colour the tree.
(((a:1,b:1[&NHX:C=2.2.2]):2[&NHX:C=2.2.2], (c:1[&NHX:C=100.100.100], d:1[&NHX:C=100.100.100]):3[&NHX:C=10.10.10]):1);

The next thing is a viewer that accepts NHX format and here I used 'Icytree' here
The alternative format a tree viewer might take is:
(((a:1,b:1[&&NHX:C=2.2.2]):2[&&NHX:C=2.2.2], (c:1[&&NHX:C=100.100.100], d:1[&&NHX:C=100.100.100]):3[&&NHX:C=10.10.10]):1);

I noted you'd used ETE3 toolkit website for your tree. Just to mention ETE3 has two levels,

the underlying raw Python library
the website GUI

The underlying raw Python is cool, except it was written to v3.6 because you combine it with your own code, the website is alright, but I think there's better.
Using NHX allows the branch width to be set as well:
(((a:1,b:1[&NHX:C=2.2.2:W=3]):2[&NHX:C=2.2.2:W=3], (c:1[&NHX:C=100.100.100:W=3], d:1[&NHX:C=100.100.100:W=3]):3[&NHX:C=10.10.10:W=3]):1);

Sadly you can't set the font (which you obviously can in ggtree
NHX normally used for auxiliary information, such as species name, or a particular phenotype.
In summary, I do prefer NHX if you are coding over R.
Note For tree colouring you have to set the use NHX format BTW in the options, a tree browser will not do it automatically.

The basic ggtree first step is to make the tree ...
library("ggtree")
tree <- read.tree(text = "(((a:1,b:1):2, (c:1, d:1):3):1);") # my bracketing might be wrong
p<-ggtree(tree)
p+geom_tiplab()

Then use ggplot2 style functions to overlay the presentation, I'm assuming you do understand ggplot2. What the code you showed was doing was individually changing every font and the position of every branch in the tree. This is definitely not cool code and I wouldn't consider it further, i.e. its not even worth bothering with, that ain't how graphs nor trees are supposed to be represented.
In any case, overall ggplot2 style code is always fiddly and to be honest I gave up using it, but you increment/or layer the changes you want as the code you supplied describes. Honestly I don't want to get into the fiddly incremental layered changes that ggplot2 graphs require.
ETE3 is good and does do coloured trees, which are described here, http://etetoolkit.org/docs/latest/tutorial/tutorial_drawing.html. The only hassle it was written for Python 3.6 - that long ago, but you can access it via conda and set the Python version.
ETE4 will be very cool but its badly delayed in its release date, they did get funding for this so the delay is not cool, maybe never who knows. Getting it up to Python 3.10 would be useful at the very least.
My personal take is ETE3 is useful for automating the tree output so you can manually assess the iterations, its quite critical for e.g. testing and debugging. Whether you are turning that into publication quality - which is universally how graphs are handled - I personally find excessive. The difference is that complex graphs represent loads of information which is not always immediately obvious to a casual observer, therefore gamma (transparency), plot type loads of stuff are relevant. A tree is a tree and 'dressing up a tree' makes not a jot of difference for an experienced observer. If someone wasn't quite clear about a tree, they'd just ask you for the tree file and file it through 'Figtree'.

See change colors or line types of arbitrarily selected branches in ggtree documentation for details. Here's my R code for applying a color to a specific branch in ggtree:
suppressPackageStartupMessages(library(ggtree))
# read tree data
tree <- read.newick('tree.nwk')
# create a dataframe of the default color (black) for branches
branchColor <- data.frame(node = 1:Nnode2(tree), colour = 'black')
# set a color for the specifc branch
branchColor[7, 2] <- 'blue'
branchColor[3, 2] <- 'orange'
branchColor[6, 2] <- 'purple'
branchColor[c(1, 2), 2] <- 'red'
# plot tree
ggtree(tree, size = 1) %<+%
  branchColor +
  aes(colour = I(colour)) +
  geom_tippoint(size = 2, color = 'black') +
  geom_tiplab(size = 5, color = 'black') +
  geom_treescale(x = 0, y = 0, width = 0.5, fontsize = 5, linesize = 1)

The tree looks like this:


text:
How to make a violin plot around quasirandom dots
I want to have a plot violin plot around quasirandom dots - like below

This is my data
> str(long1)
Classes ‘data.table’ and 'data.frame':  116 obs. of  4 variables:
 $ gene    : chr  "WT1" "KRT7" "MUC16" "EPCAM" ...
 $ variable: chr  "1016" "1016" "1016" "1016" ...
 $ value   : num  5.82 8.41 11.22 7.21 5.4 ...
 $ group   : chr  "Tumour" "Tumour" "Tumour" "Tumour" ...
 - attr(*, ".internal.selfref")=<externalptr>
>

By this code I get plot which does not have any violin
library(ggplot2)
library(ggpubr)
library(ggbeeswarm)
p <- ggplot(long1, aes(x=group, y=value,fill=group))
p +geom_quasirandom(method = 'quasirandom',alpha = 0.2, width = 0.2,aes(colour = variable))+ labs(title="Supplementary Figure 6B", y = "Expression(logcounts)")+theme_classic()+
  theme(legend.position="right") +
  theme(legend.title=element_blank())+facet_wrap(~gene, scales = "free")+theme(strip.background = element_blank(),strip.text.x = element_text(size = 12,face="italic"))

I also tried adding geom_violin like
p <- ggplot(long1, aes(x=group, y=value,fill=group))
p +geom_violin()+geom_quasirandom(method = 'quasirandom',alpha = 0.2, width = 0.2,aes(colour = variable))+ labs(title="Supplementary Figure 6B", y = "Expression(logcounts)")+theme_classic()+
  theme(legend.position="right") +
  theme(legend.title=element_blank())+facet_wrap(~gene, scales = "free")+theme(strip.background = element_blank(),strip.text.x = element_text(size = 12,face="italic"))

But gives a weird plot not similar to what I want

Any help please
EDITED
Thanks to @haci
By this code I get this plot but as you see the point color inside the violin is not really similar to what I want
Any suggestion to get a similar color scheme for dots inside the violins
p <- ggplot(long1, aes(x=group, y=value,fill=group))+guides(fill = "none")
p +geom_violin(alpha = 0.1)+scale_fill_manual(values=c("#fcf4e0", "#fcf4e0"))+geom_quasirandom(method = 'quasirandom',alpha = 0.2, width = 0.2,size=2,stroke=0.5,aes(colour = variable))+scale_color_hue(l = 5, c = 150)+ labs(title="Supplementary Figure 6B", y = "Expression(logcounts)")+theme_classic()+
  theme(legend.position="right",legend.text=element_text(size=12)) +
  theme(legend.title=element_blank(),legend.text.align = 0,text = element_text(family = "Arial"),
        axis.title.x = element_text(size = 0.00001),
        axis.title.y = element_text(size = 17))+theme(axis.text = element_text(
          color="black",
          size=12,
          face=1)
        )+facet_wrap(~gene, scales = "free")+theme(strip.background = element_blank(),strip.text.x = element_text(size = 17,face="italic"))+ theme(axis.line = element_line(colour = 'black', size = 0.8))

You are almost there, if you set the alpha parameter within your geom_violin() call (to sth like 0.1) you should be able to get a plot similar to the example plot.
The fact that you only have two data points for the Stroma group and these two points having similar values, you will end up getting a horizontal-line like violin for this group. Moreover, you do not have that many data points in general anyway so you might consider not using violins at all, especially if the outliers in Tumour group is coming from the same (maybe somewhat problematic?) sample, that might need to be dropped from the analysis. The reader should be able to grasp the essence of the data with our without violins (violins could help making the outliers obvious though).

Thanks @haci . In edited version of my post I am struggling with the same color scheme, please do you have any suggestion here?

What you mean @M_

@Exhausted Don't edit your question and branch off, ask a new question as a follow up (by mentioning this question) instead. When you edit questions that have been answered, you confuse future users by removing context. Scientific forums are not your personal help desk.

@Exhausted, you can try increasing the alpha parameter in your geom_quasirandom() call (or remove that parameter altogether).


text:
Can we do NGS library preparation using UMI with large amount of DNA input?
We know that in next-generation sequencing (NGS), the unique molecular identifier (UMI) can reduce or eliminate sequencing or PCR errors and result in very high accurate data. Therefore, UMI is widely used in rare mutation detection. A typical application is cancer early detection or minimal residual disease (MRD) detection that try to find the very rare DNA from cancer cells in the blood.
However, we also know that the cell-free DNA in the blood is very low, typically 5 ng in 1 mL blood. And it usually takes about 10 mL blood to get about 50 ng DNA to get sequenced with UMI. Limited DNA input leads to many false negatives.
So I came up with an idea: What if I got a sample with plenty of DNA? Can I detect rare mutations using UMI in this sample?
The UMI library preparation protocols I have met all use a small amount of DNA input, say 50ng. If I have 1μg DNA, do I have to cut it to 50ng to prepare the UMI library?

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

Do you mean that the variant allele frequencies (VAFs) of the variants will be lower than 1% for the very rare DNA from cancer cells in the blood?

There are around 3000 copies of haploid genome present in 10 ng of DNA, approximatively (50-)60 ng of cell free DNA will be required for a sensitivity of 0.01% (one rare event in 10,000 molecules),  moreover,  in blood, average cfDNA concentration in blood of cancer patients, varying between 0 and 1000 ng/mL, with an average of 180 ng/mL

The question is

"Do I add 20-fold more DNA than the kit states?"

Its not my area - this is wet-lab - but there will be a tolerance however 20-fold is excessive (>1 log), i.e. no it doesn't help and could cause (expensive) downstream problems. The people to contact are the support for the given kit, there are unlikely to support your idea.
The answer isn't wet-lab but almost certainly is bioinformatics and detecting low frequency somatic driver mutations. The Broad Institute worked on this problem for years, its nothing new, i.e. identifying driver mutations in highly mixed sample comprising healthy somatic tissue and cancer cells.
I've tracked their work and have a reasonable understanding of what they do and why: oncology isn't my field. Your in-house bioinformatician should know this stuff.
One suggest is using MuTect from the Broad Institute, I assume there have been many custom solutions since then. Their famous product is GATK for VCF, I don't think thats the solution here.
Again cancer isn't my area BTW this is just common knowledge within bioinformatics, i.e. knowing what someone else is doing and why is part of the knowledge base.
Please do try to retain a bioinformatics focus for future questions.

Whole genome sequencing uses between 50 ng and 1 μg of high-quality DNA (A260/A280: 1.7–2.0) for PCR-amplified or non-amplified libraries if your focus is precision oncology as described in the following open access review in Seminars in Cancer Biology paper by Meggendorfer et al 2022

Do you mean the following type of technology?

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


text:
Can FASTA data be submitted directly into eCTD submission to the FDA?
We are getting sequence data from microbiota samples as part of a FDA submission identifying changes in Micobiota from skin samples over time and compared with different treatments.  the data will be first presented in FASTA formate then analyzed for OTUs and differences.  The submission has to follow eCTD format.

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:
Problem with sequence alignment
I'm trying to find if the peptides in column A are present within the longer peptides of column B.
Pairwise alignment tools take the data from both columns as one string (i.e., TIIDYTNNHYTNNHLEII... and PFYMNITGYTDAYNTDIIDYIIDYTNNHLEIITLFN...) which leads to a number of peptides from column A left unaligned even though the peptide is present within one of the longer peptides in column B.
I would like to find a method that works for detecting such peptides.

Easy, few details are you stuck with a CSV format or prepared to copy and paste column A to one file and column B to another. Its solvable either way.

I don't think that you want an alignment so much as you want to test for the presence of an exact substring (column A) within larger strings (column B).
There are simple excel utilities that do this, for example the SEARCH() function. There is no need to use alignments.


text:
Basic RNA Differential Expression in R
I have two matrices, one for individuals before treatment and one for the same individuals after treatment. Both matrices are raw read counts of RNA expression.
        Treated_1  Treated_2  Treated_3  Treated_4  Treated_5
RNA_1   105        283        64         155        51
RNA_2   359        27         47         348        84
RNA_3   99         10         89         345        77
RNA_4   48         100        77         74         83

        Untreated_1  Untreated_2  Untreated_3  Untreated_4  Untreated_5
RNA_1   100          200          50           130          33
RNA_2   200          10           30           300          75
RNA_3   90           10           60           320          60
RNA_4   66           50           33           60           55

How can I tell which RNA's are significantly differentially expressed in R using edgeR or DESeq2? I'm thinking first normalize using TMM then compare the mean of the normalized values for treated vs. untreated?

See here for a Differential Expression guide which discusses carrying out differential expression using DESeq2:
https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
It starts with a standard "quick start" workflow, and goes into quite a lot of detail explaining all the different things that can be done:
dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design= ~ batch + condition)
dds <- DESeq(dds)
resultsNames(dds) # lists the coefficients
res <- results(dds, name="condition_trt_vs_untrt")

I'm thinking first normalize using TMM then compare the mean of the normalized values for treated vs. untreated?

That is a bit too simple. You need to analyse your experiment as a paired comparison. You need to combine the treated and untreated count matrices into one and form factors for individual and treatment:
Counts <- cbind(Counts.treated, Counts.untreated)
Individual <- gl(5,1,length=10)
Treatment <- gl(2,5,length=10,labels=c("Treated","Untreated"))
Treatment <- relevel(Treatment, ref="Untreated")
design <- model.matrix(~Individual+Treatment)

Then follow a standard workflow:
https://bioconductor.org/packages/release/workflows/vignettes/RnaSeqGeneEdgeRQL/inst/doc/edgeRQL.html
library(edgeR)
y <- DGEList(counts=Counts)
y <- calcNormFactors(y)
y <- estimateDisp(y,design)
fit <- glmQLFit(y,design)
test <- glmQLFTest(fit)
topTags(test)


text:
Retrieve chromosome number with just RSID?
I have ~500,000 rsids - but want which chromosome each SNP resides on. How can I do this?
MarkerName A1 A2
rs2326918 a g
rs7929618 c g
rs66941928 t c

I want the chromosome number associated with each MarkerName. So the output would be like this:
MarkerName A1 A2 Chr
rs2326918 a g     6
rs7929618 c g     6
rs66941928 t c    6

How can i do this?

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

Thanks for the clarification @workhorse

I'd recommend consulting the RefSNP documentation https://www.ncbi.nlm.nih.gov/snp/docs/RefSNP_about/

@Workhorse thanks for your question. Please consider upvoting/accepting the answer by NiekdeKlein

Download dbsnp VCF file from https://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/. This contains rs IDs and chromosome info (among other info). Can match your rs IDs to the rs IDs in dbsnp VCF, and extract chromosome from there.
See https://www.ncbi.nlm.nih.gov/variation/docs/human_variation_vcf/ for more information on the dbsnp files.


text:
Aligning scRNA-seq fastq to .bam without cell barcodes
I would like to convert the fastq files from this dataset: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE98679 to .bam files.  The group did not provide the scripts they used to to their alignment, I know they already posted their reads but I need the .bam files for my pipeline.  Is it possible to do this without knowing the cell barcodes?  If so with what tools?

The linked dataset is based on the Smart-seq2 technology. This is a plate-based assay in which individual cells are first sorted into individual wells of a microwell plate via FACS. The cells are then individually lysed and libraries are prepared. Since this happens per well one assigns a single sequencing barcode to each well/cell, hence no further cellular barcoding as in droplet-based methods such as 10x Chromatium is required. Long story short, each fastq file that one obtains from the sequencing represents one individual cell. Hence, one can use any standard RNA-seq alignment or quantification technique such as STAR or salmon to obtain alignments/quantifications. One would then need to parse the individual quantifications into a single count matrix. If one used aligners such as STAR this could be done with high-performance quantifiers such as featureCounts. If using pseudo- or selective aligners such as salmon or kallisto I would use tximport to aggregate the per-transcript quantifications to the gene level while importing into R to get a gene-level count matrix. As each file is a single cell there is no need for more complicated solutions such as CellRanger which do not apply here.

@foobar.2   What other sequencing technologies is this true for?  What about InDrop for instance?


text:
How can I remove FASTA sequences with an asterisk inside the sequence?
I have several protein sequences with an asterisk in the middle of the sequence, and I'd like to remove those sequences.
In the example below, I have three sequences:

From this example, I would like to obtain only the following output:

Why do you want to remove sequences with * in them?

Try seqkit grep
seqkit grep --by-seq --invert-match --pattern '*' seq.fa

--by-seq search by seq
It will output sequences without an asterisk.

This technically answers the question, if we assume OP is not in an XY problem and removing valid sequences that they don't understand (which is highly probable given that they call them "sequences with an asterisk" and not "sequenced containing stop codons".

It’s  good you brought up the question.
* are stop codons
Don't simply parse out the * symbol. If done incorrectly, it could constitute 'tampering' and that’s a red line in biology.
The * means a stop codon in the middle of a protein sequence. There are quite a number of reasons why there is an *, i.e., a stop codon in the sequence, for example, in no particular order:

the open reading frame is out of frame which is part of the natural molecular biology, i.e., there can be multiple reading frames in a gene
the open reading frame is out of frame due to a sequencing error and poor annotation (which is common);
the wrong reading frame was used (very, very common);
the locus isn't translated, e.g., it’s an intron or UTR (untranslated region)

The list goes on, e.g., it could be expression regulation ...

If it is not an intron/UTR the way you might 'remove' * is by changing the reading frame, e.g., try using the second codon position (then the third) as the opening reading frame: if the * 'disappear' it was very likely just a reading frame issue. How to do that is a separate question. The resulting amino acid sequence will of course be very different.

There's a second possibility if the reading frame is correct, which is to treat each locus between the * as a separate protein. That way the * is 'removed' because their a multiple proteins within the locus. In that case, a \n would replace the * (Mac and Linux), i.e., carriage return.

I don't think they're looking to tamper, but they probably don't understand the significance of the *. Your point should help them avoid an XY problem if they're in one - I hope they read your answer before proceeding in whatever direction they're headed.

You could use Biopython for this. It makes FASTA parsing really easy.
'''
Python version: 3.10.0
Biopython version: 1.79
'''
from Bio import SeqIO
for seq_record in SeqIO.parse("input.fasta", "fasta"):
    if "*" not in (seq_record.seq):
        print(">",seq_record.id)
        print(seq_record.seq)


text:
System specifications for NGS data analysis
I have a 3.7 GB whole genome data of a eukaryote, for which genome assembly, gene prediction, and annotation steps have to be performed. In some time I would also need to analyze the transcriptome data of the organism. Can someone please let me know what are the system requirements for carrying out the aforementioned NGS data analysis? 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.

"3.7 GB whole genome data" what kind of data. sequencing? what platform? what tech specs?


text:
Good molecular docking software?
I'm looking for some good molecular docking software (preferably, one that comes with a Python API). I obviously found Vina and that's been quite helpful but I was wondering if there's anything out there that's more accurate and reliable ? Since there's so much software and publications out there, I thought someone here might be aware of some things that may be overlooked.
Thanks in advance !

The Schrödinger suite has among the nice tools, Glide, a docking tool, which generally outperforms other docking tools in benchmarks, but it is expensive. Generally however there will be a couple of licences that could be borrowed floating around in a given department —always worth asking around. However, docking is inaccurate —full stop— thus is a numbers game, so if say one in ten is legit, then you buy 50 compounds for testing.

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

The Bot has a good point: it's a topic of ample discussion heavily benchmarked (occasionally with biased benchmarks): nobody wants to write a review here. However, focusing on a particular task can help greatly get an answer on how to improve your docking task as good experimental setup can change things a lot and will beat Autodock Vina under a generic setup: Say you want to dock a library in the adenine site of a kinase vs. you want place a GTP instead of a non-hydrolysable analogue vs. you want to dock the native ligand in an enzyme.

Thanks very much for all the replies. I come from a wet-lab based background so I don't like docking myself but some reason it appears to be a requirement for workflow validation. So I thought if I'm going to add something I don't have confidence in, I might as well add the most accurate one I can.
The project is meant to be easily to modify to adapt to new scenarios so that's why the generic question.

"Method A is better than B" is always tricky when they are utterly different and will always be rather biased (by both the experimenter and reviewer 2) and one could argue that if method B were done differently taking into account X it would win.
The major improvement with docking one can get is feeding it more information, which if the method A uses already then that argument is rather strong.

Examples of improvements making the box restrictive to the active/binding/hinge site (there are too many question here on the lines of "Vina is telling me 30 Å width is too big but I am using the whole xstal cell").
Constraining based on known binders.
Dealing with isomers (stereo- and tauto-).
Correcting the protonation of residues and of the ligand.
Etc.

However, here it sounds like you want to do docking just because other have and it does not enrich your story —then don't muddle it: empirical results (xstal, assays etc.) trumps docking hands down. Those papers may just have the data handy as they did it in the planning step even if the docking is shown as a latter step in the final paper...


text:
Common tools for generating results and visualizations for bioinformatics and genomics publications
I am a computational scientist with background in computer science and engineering however very new to the field of bioinformatics and computational genomics. Looking at papers in nature, PLOS, Cell, PNAS etc. within the topic of computational genomics, I see a very similar set of very sophisticated graphical results, ranging from clustered/sorted heat map representation of gene expressions to log2 fold changes, volcano plots, circular graph network visualizations etc. Coming from an engineering background I would often generate customized tools for visualizing results myself as part of each paper, albeit not at the sophistication level of results I see in bioinformatics. But here, it seems there is a set of tools standardized where many are using them to display their results. Is there a summary of these result generation/visualization tools along with the references to their implementation/usage (opensource I assume)? Something like that could be very helpful to newcomers like myself. Thanks in advance.
Below is an example /elaboration based to guide the question further into specifics:
Consider the example paper:
https://doi.org/10.1016/j.cell.2018.04.018
This paper contains a principle set of results/visualizations I see almost universally in many nature/cell computational genomics papers. I am wondering if there exists an integrated library either in python, R (or Matlab for example when they create toolboxes and integrate a commonly used subset of tools in an area) that everyone taps into for generating these type of results. The "principal set" that I am referring to include:

Genome browser track snippets: Figure 1A
log2Fold x-y scatter plots : Figure 1B
Hierarchically clustered/sorted gene exp heatmaps: Figure 1C,1F
Venn diagrams : Figure 4F
Hi-C type square ligation heat-maps: Figure 6A

Figure 1A: These genome browser snippets are either plain screenshots from applications like IGV or created from within R with software like gviz.
Figure 1B: This looks like base R plotting routine. Like plot(x,y) and then some coloring and legending. No particular "tool" here other than really just base R or ggplot2.
Figure 1C/F: Heatmap packages include ComplexHeatmap and pheatmap. I personally use the former as it is versatile with extensive documentation.
Figure 4F: Venn diagrams, just google "Venn" and "R" for several packages that do it in R. Same for other languages. There are plenty of packages. No preference.
Figure 6A: These plots are often returned by specialized HiC software such as HiCexplorer which is powerful with lots of documentation.


text:
Statistical significance of pairwise protein sequence alignment
I have identified a conserved region in a protein I'm studying, by doing multiple sequence alignments against a dozen of its orthologs. I've then found that this region also aligns relatively well with a functional region (binding motif) in a completely unrelated protein, and I'm trying to estimate whether it's reasonable to believe the two regions might be bound by the same protein. blastp doesn't identify this match, would I then assume that this is insignificant?
Let's say I have a 15 aa long peptide sequence X, which I align against a longer 200 aa protein sequence Y. I find some best match (gaps allowed).
What statistical test could I now perform to answer the question:
Is the match I found for sequence X in sequence Y significantly better than what I would expect to find between X and a random 200 aa sequence Z?

Just formalise the answer, its a straight use of the Blast E-value and I wouldn't personally advise against anything beyond that (explained below) within molecular evolution.
So the Google gives the definition

The BLAST E-value is the number of expected hits of similar quality
(score) that could be found just by chance. E-value of 10 means that
up to 10 hits can be expected to be found just by chance, given the
same size of a random database.

Thus the smaller the E-value the better the hit.
Using a local blast server will provide a specific query verse reference E-value.
blastp -query gene.fa -subject genome.fa

The broader question concerns the whether the Blast hit is a genuine homology and that is very complicated. This is simply because amino acid loci could have arisen by convergent evolution rather than direct descendancy which are homologues and this can become controversial.
When trying to assign a function to a amino acid locus this starts becoming relevant because e.g. the selection pressures of convergence could very difference.
Identifying homologues would require a phylogenetic tree which is compatible with the established phylogenetic tree, e.g. for house-keeping genes as a starting point. It gets complex quickly beyond that. Functional domains within the Blast variant (Blast-Psi) could be of use, its field specific however.
Thus the simple answer is to use Blast E-value (bit-score) and it is a universally recognised metric, the reader can then make up their own mind what the value represents. You could simply simple state the definition of the E-value in the report.

Thanks for the recommendations, this is very clear. I was able to tweak the blastp parameters and got an expect value of 0.15, but like you point out I'd have to make up my own mind whether this is relevant or not (a least it's below 1 ;)). Also, in this specific case I'd expect the two regions to have arisen independently via convergent evolution, so no shared phylogeny expected.

Good its sorted, thanks


text:
Determination of transmembrane domain limits using multiple prediction tools
Integral membrane proteins are proteins that have at least one transmembrane (TM) segment.
I performed predictions of TM domains using three tools: ProtScale, ΔG predictor 1.0, TOPCONS.
For example, the results in terms of aminoacid positions for the TM1 domain were as follows:

ProtScale: 37-55
ΔG predictor: 37-55
TOPCONS: 36-56

The limits of the predicted TM domains can be used for drawing the topology of this protein.
For this purpose, we need to have representative limits from these results.
Statistically speaking:

which values should be chosen?
should we calculate a mean : (37+37+36)/3=36.67 knowing that these values correspond to aminoacid positions and 36.67 doesn't mean anything?

How can I use these results (37, 37, 36) to determine the limit corresponding to the 1st aminoacid of TM1 segment?

How big is the protein? Normally TMDs are at the 3' end of the protein.

The protein has 426 aminoacids. TM1 is the first TM segment, from the N-terminal domain.

The variation you are seeing are minor differences in estimates in the transmembrane domain, these are usually modelled by hydrophobicity and are accurate predictions. If you look at the secondary/tertiary structure you might see alpha helices which will 'burrow' into the membrane. The region in this particular question is a small domain however, but the position is 'about right'.
Thus question is whether 1 amino acid at 5' and 3' site is inside or outside the membrane - I think this is overthinking. Edge-effects are common in any estimate and 1 amino isn't bad variation at all. To double check have a look at the Kyte-Doolittle index on the Expasy website (which is one of the first), and TMHMM. Finally, you need to check the crystal structure if available.
What does puzzle me is the apparent absence of a 3' GPI anchor, but I dunno the biological system so ain't for me to comment.


text:
Converting Chip-seq fastq format ({$x_k$} list of seq segments) to genome browser {sequence $x_i$, count-number $y_i$} format
Very new to bioinformatics. I am trying to figure out the steps required to convert a chip-seq fastq file into a genome browser track, with the end goal of being able to access the sequence (x) vs. activity level/count (y) data associated.
Any genome browser activity/count track data has the following functional format

$y=F(x)$
Here x is the sequence axis, and y is the activity/count axis. Typically the x(i) is the nucleotide associated with the precise genomic coordinate i (an integer number) depending on the genome used to align (human, mouse, etc.). What I need is a format allowing direct access to the nucleotide representation x(i), associated with the index i, as well as the count levels at y(x(i)). What is the closest data file type that I need? And how do I get to it?
In other words, what are the steps for a fastq file associated with a chip-seq data to be transformed to this x(i),y(i) format, or best yet if there is an example Python or R code out there that could illustrate the steps of going from file1.fastq to something like to two files file1_x_sequences and file2_y_counts? I only know the first step is an alignment to a genome, but stuck on the next steps to complete the picture. Thank you.
Ultimately I need a list of sequences (tags) lets say 250bp long, and a count number that represents the chip-seq enrichment for the part of the genome (the 250bp region). Further reading, indicates I need to detect/call "peaks" with MACS2/homer to generate a BAM file and ultimately a Fasta file for the sequence lists, but I don't know yet how the 'counts' come in (in what file format and where from). Any hints there would complete the picture I think. Ultimately this is training data for x=250bp sequence to y=count/chip-seq enrichment activity model.

I have no idea about what you are actually trying to achieve, but getting these browser tracks is simple. First you align your fastq files to a reference genome, e.g. with bowtie2, this returns a SAM file. Sort the SAM file, convert to BAM and sort + index it with samtools. Then use something like bamCoverage from the deeptools suite to get your browser tracks.
A word of caution. Since you seem to be new to OMICS and bioinformatics in general: NGS data, especially applications that are as noisy as ChIP-seq have lots of pitfalls. Be sure to go through the literature to find references for what you are trying to develop, at least the basic idea and discuss with experienced folks whether it makes sense. It sounds like you want counts of the peak summits, is that correct? In any case, be sure not to reinvent the wheel and to do a meaningful analysis. Discussion with experienced peers is key, that may save you from pitfalls and rookie mistakes. This entire choice of words "sequence axis", "activity level" is very uncommon, I have never heard of this before despite working on NGS data for years, for what it's worth.


text:
How to create Seuratobject using 3 samples
I have total of 6 scRNA seq data samples which can be divided into group A and B. I would like to create Seuratobject A and B by combining 3 samples counts data for each group. Eventually I would like to see UMAP with colored by group A and B. I tried many ways but failed. I can see UMAP with 6 colors representing 6 samples. Could you please help me with this issue?
Please find below my code for your reference-
A1.data <-Read10X(data.dir = "C:/Users/mad1188/Box/SequencingData/cellranger_count/-A-7-1/filtered_feature_bc_matrix")
A2.data <-Read10X(data.dir = "C:/Users/mad1188/Box/Sequencing Data/cellranger_count/-A-7-2/filtered_feature_bc_matrix")
A3.data <-Read10X(data.dir = "C:/Users/mad1188/Box/Sequencing Data/cellranger_count/-A-7-3/filtered_feature_bc_matrix")
A1 <- CreateSeuratObject(counts = A1.data,
project = "A7")
A2 <- CreateSeuratObject(counts = A2.data,project = "A7")
A3 <- CreateSeuratObject(counts = A3.data,project = "A7")
head(x = colnames(x = A1))
A1 <- RenameCells(object = A1, add.cell.id = "A7")
head(x = colnames(x = A1))
head(x = colnames(x = A2))
A2 <- RenameCells(object = A2, add.cell.id = "A7")
head(x = colnames(x = A7))
head(x = colnames(x = A3))
A3 <- RenameCells(object = A3, add.cell.id = "A7")
head(x = colnames(x = A3))
A1[["percent.mt"]] <- PercentageFeatureSet(object = A1, pattern = "^MT-")
A2[["percent.mt"]] <- PercentageFeatureSet(object = A2, pattern = "^MT-")
A3[["percent.mt"]] <- PercentageFeatureSet(object = A3, pattern = "^MT-")
p1 <- VlnPlot(object = A1, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
p2 <- VlnPlot(object = A2, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
p3 <- VlnPlot(object = A3, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
CombinePlots(plots = list(p1, p2, p3), ncol =1)
p4 <- FeatureScatter(object = A1, feature1 = "nCount_RNA", feature2 = "percent.mt")
p5 <- FeatureScatter(object = A1, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
p6 <- FeatureScatter(object = A2, feature1 = "nCount_RNA", feature2 = "percent.mt")
p7 <- FeatureScatter(object = A2, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
p8 <- FeatureScatter(object = A3, feature1 = "nCount_RNA", feature2 = "percent.mt")
p9 <- FeatureScatter(object = A3, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(p4, p5, p6, p7, p8, p9), ncol =2)
A1 <- subset(x = A1, subset = nFeature_RNA > 200 & nFeature_RNA < 8000 & percent.mt < 20)
A2 <- subset(x = A2, subset = nFeature_RNA > 200 & nFeature_RNA < 8000 & percent.mt < 20)
A3 <- subset(x = A3, subset = nFeature_RNA > 200 & nFeature_RNA < 8000 & percent.mt < 20)
p10 <- VlnPlot(object = A1, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
p11 <- VlnPlot(object = A2, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
p12 <- VlnPlot(object = A3, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
CombinePlots(plots = list(p10, p11, p12), ncol =1)
p13 <- FeatureScatter(object = A1, feature1 = "nCount_RNA", feature2 = "percent.mt")
p14 <- FeatureScatter(object = A1, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
p15 <- FeatureScatter(object = A2, feature1 = "nCount_RNA", feature2 = "percent.mt")
p16 <- FeatureScatter(object = A2, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
p17 <- FeatureScatter(object = A3, feature1 = "nCount_RNA", feature2 = "percent.mt")
p18 <- FeatureScatter(object = A3, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(p13, p14, p15, p16, p17, p18), ncol =2)


text:
plotting gene expression after EdgeR DE analysis using RUVg (RUVseq) covariates
I have used the empirical RUVg method (from RUVseq) to estimate the unwanted variation of my dataset (consisting of several public datasets analysed together, with controls and case samples in different datasets.
Afterwards, I used EdgeR to perform my DE analysis with 2 RUV covariates:
designRUV <- model.matrix(~0 + condition4 + W_1 + W_2, data = pData(set2))

y <- DGEList(counts=counts(set2), group=pheno$condition4)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, designRUV)
y <- estimateGLMTagwiseDisp(y, designRUV)
fit <- glmFit(y, designRUV)

I got my results after making all the contrast that I needed but I would like to plot the expression of some genes of interest. One plot per gene with all the conditions in condition4 has 9 categories (expression level or norm counts on the Y axis).
My doubt is which counts to use. I have the normalised ones from RUVg, but it states they should only use them for exploratory analysis (maybe this is the case??). These counts as I understand are the residuals of the model, so not really sure what to do.
Another option is to use the normalised counts produced by EdgeR, but these ones doesn't account for the unwanted variation (W_1, W_2 covariates) (right???)
Any advice or nice solution?? :=)
Thanks in advance!!

You could correct your logCPMs (edgeR::cpm(y, log=TRUE)) with the covariates from RUVseq using limma and then use these for plotting:
set2 <- RUVg(...)

designRUV <- model.matrix(~0 + condition4 + W_1 + W_2, data = pData(set2))

y <- DGEList(counts=counts(set2), group=pheno$condition4)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, designRUV)
y <- estimateGLMTagwiseDisp(y, designRUV)
fit <- glmFit(y, designRUV)

# correction with limma
logcpm <- cpm(y, log=TRUE)
design <- model.matrix(~0+condition4, pData(set2))
corrected <- removeBatchEffect(logcpm, design=design, covariates=as.matrix(pData(set2)[,c(W_1,"W_2"])

Thank you @alderaan_shot_first. Can I ask why you used a "-" in the model.matrix? Also, after reading your answer I have been checking the removeBatchEffect documentation and when you put the covariates argument you are setting all the phenoData of set2... Don't would be covariates=as.matrix(pData(set2)[,c(W_1,"W_2"]) ?. Thank you!!!

was a typo, sorry. As for the matrix, yes it is fine with the code you suggest. feel free to make an edit.


text:
How to extract DNA sequence from browser track input files (BigWig, bed etc.) files
A BigWig file for something like Chip-seq, or ATAC data has the following format for each line i:
chrom_i - start_i - end_i - value_i

I need to convert for each line the (start_i -> end_i) reference to the actual DNA sequence 'ATATCCTAG...' so that at the end I have a list of DNA sequences like 'ATATCCTAG...' and the Chip-seq activity value associated with that sequence segment. What are the ways to do this? For example, what tool(s) can I use to generate the (sequence, value) pairs from browser tracks for the Chip-seq or ATAC etc. And in fact, going from BigWig browser track files for Chip-seq the most direct way to produce this kind of (seq,value) data? Thank you.

[Here is the link.][1]. Instructions are very clear at the website
[1]: https://bedtools.readthedocs.io/en/latest/content/tools/getfasta.html

Ah, perfect @Supertech! Thank  you, both answers helpful... the link included here is the exact answer I was looking for: bedtools getfasta converts index range start_index-end_index into it's equivalent sequence of nucleotides indicated in the input *.fa file. I just need to figure out how to create the "fa" file for an entire genome like mm10, where it's usually in chromosomal chunks (i.e. Chr2.fa).

I believe each chromosome is independently numbered. If you want to generate a single genome file with all the chromosomes, download individual chromosomes and simply concatenate them. For example using Linux command "cat". >cat chr1.fa chr2.fa ...

It looks like they may have single file. I have not downloaded to check it though. It has been long that I did not work with mouse genome. http://hgdownload.cse.ucsc.edu/goldenpath/mm10/bigZips/

Yes, I am trying to think why/how that would work: I guess the *.fa file has the chrX indicator at each line which would be matched along with the index (I assume a complete index needs to include chromosome number as well)! Thank you again.

Try samtools faidx. Check bedtools manual. https://bedtools.readthedocs.io/en/latest/content/general-usage.html

Bedtools can take a list of coordinates and a fasta and return a multifasta of the specified regions.


text:
How to identify which isoforms of a gene are actually expressed in my data?
I am new to bioinformatics and have been assigned the task of discovering which isoforms of a certain protein-coding gene are present in the RNA-seq data that I have. There are several isoforms of this gene listed on Entrez & Ensembl and I want to know which of these isoforms are expressed in the RNA-seq data that I have been given.
I am not sure how to begin to tackle this issue so If anyone has any ideas or advice, it would be greatly appreciated.

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

Hi @Sawason, I've provided an answer to your question, but I'm not totally sure if it's what you're looking for. It might help if you edit your question to provide some additional information, like what species you're working in, or some accession IDs if you're working with publicly available data, for example. It's tough to give specific advice for a question that's this general when we don't really know what you're trying to do.

Tools like kallisto or salmon can perform a process call pseudo-alignment that matches reads from RNA-seq data to transcript isoforms present in a reference genome.
You can then look at the output of these tools to see which isoforms for a given gene appear to be present in your RNA-seq data.
These tools answer the question "how abundant is each transcript in my data?".
But if you want to answer a question like "are there differences in which transcripts are present between experimental conditions?", I'd look for tutorials on differential expression analysis, differential splicing analysis, or differential transcript usage.

If you do as @james-hawley rightly suggests, keep in mind that pseudo-alignment could detect  a lot of transcripts with low TPM that are actually false positives, as reported here (it's a paper focused on scRNA-seq, but many statements are also valid for bulk). So, you may want to try to remove false positive transcripts, that will likely result with very low read count and TPM. Salmon's decoy function could help here.
As an alternative to pseudo-alignment, you could run StringTie on your BAM file to detect transcript isoforms, using a reasonable threshold of coverage over splice junctions.