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:
scRNASeq expression matrix with decimal values
I am trying to replicate some results of a scRNASeq experiment and, when I looked at the data provided by the author, I noticed that some of the counts in the expression matrix are represented as decimals.
Would this be because the authors uploaded a normalised count matrix or for other reasons? Also, please note that the authors used the 10X genomics cellranger pipeline for the processing of raw data.
I have provided a brief snippet of the data frame here:
structure(list(P1TLH_AAACCTGAGCAGCCTC_1 = c(0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), P1TLH_AAACCTGTCCTCATTA_1 = c(0,
0, 0.314759545319035, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0.314759545319035, 0, 0), P1TLH_AAACCTGTCTAAGCCA_1 = c(0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1.18024827050334, 0, 0, 0, 0, 0, 1.18024827050334,
0, 0, 0), P1TLH_AAACGGGAGTAGGCCA_1 = c(0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0.644905372424165, 0, 0, 0), P1TLH_AAACGGGGTTCGGGCT_1 = c(0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.780522072779651,
0, 0, 0), P1TLH_AAAGCAACAGTAAGAT_1 = c(0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)), row.names = c("RP11-34P13.7",
"FO538757.2", "AP006222.2", "RP4-669L17.10", "RP5-857K21.4",
"RP11-206L10.9", "LINC00115", "FAM41C", "RP11-54O7.1", "RP11-54O7.3",
"SAMD11", "NOC2L", "KLHL17", "PLEKHN1", "RP11-54O7.17", "HES4",
"ISG15", "AGRN", "C1orf159", "RP11-465B22.8"), class = "data.frame")

Also, I am quite confused of the structure they used - as I don't have much experience with single cell datasets. Just to give some further context:
The colnames of the "full" dataframe are coded from P1TLH_sequence_1 to P5TLH_sequence_1. There are 5 patients in the study.
I'm assuming that each of the colnames represents a cell for a particular patient while the rows are the genes.
If this assumption is correct, how would I do the conversion from colnames in the dataframe to the "cell names"? And lastly, why are there decimals in this expression matrix?
Thank you all for your time,

It looks like you have normalized counts. My guess is that they are log-normalized with some scale factor, probably 10000.
With the whole matrix you should be able to work out the scale factor for each cell. See https://www.nxn.se/valent/2018/10/25/unscaling-scaled-counts-in-scrna-seq-data
They typical format for scRNA-seq data is to have cells in the columns and genes in the rows. So here your rownames are the gene names and colnames are the cell names.

Thanks for the response ! Do you have any idea of a DB of some sorts that is used to convert the colnames from sequence format to cellID format?

What's the cellID format?

It is not a particular format called cellID per ce. I just mean conversion from the sequences to a naming system that would be more suitable for downstream analysis / clustering.

Ok in that case any unique identifier is fine and I would leave the cell names as they are. The only thing you might like to do is extract the patient ID information that appears to be encoded in the cell name. If you use Seurat, this will be done for you when creating the Seurat object.

source: Stack Exchange

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

sort: votes


text:
Running htseq-count over BAM files
I am trying to derive an expression matrix from BAM files using htseq-count. These are bulk RNASeq BAM's by the way.
I have read the htseq-count documentation as well as samtools and figured that the following command should work:
samtools sort -on file.bam | samtools view -h - | htseq-count --type=exon --idattr=gene_id --additional-attr=gene_name - Homo_sapiens.GRCh38.77.gtf > counts.txt

The .gtf file here is from ensembl.
I have tried many variations of the above command but have failed to generate a count matrix. The error message states:
"0 SAM alignments processed"
Here are the first few lines of my .gtf file:
#!genome-build GRCh37.p13
#!genome-version GRCh37
#!genome-date 2009-02
#!genome-build-accession NCBI:GCA_000001405.14
#!genebuild-last-updated 2013-09
1       pseudogene      gene    11869   14412   .       +       .       gene_id "ENSG00000223972"; gene_name "DDX11L1";
gene_source "ensembl_havana"; gene_biotype "pseudogene";
1       processed_transcript    transcript      11869   14409   .       +       .       gene_id "ENSG00000223972"; trans
cript_id "ENST00000456328"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_nam
e "DDX11L1-002"; transcript_source "havana";
1       processed_transcript    exon    11869   12227   .       +       .       gene_id "ENSG00000223972"; transcript_id
 "ENST00000456328"; exon_number "1"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; trans
cript_name "DDX11L1-002"; transcript_source "havana"; exon_id "ENSE00002234944";

And all the genes in the text file have an expression value of 0. What exactly is going wrong in this step? Is the conversion to a SAM file the problem? I'm completely lost so any help would me much appreciated.
Thanks!
EDIT:: For anybody who sees this post cause they have a similar problem. Just note that the newer versions of htseq-count don't require sorted .bam files and, so following the editing of the .gtf file, all I needed to do was convert it to .sam using samtools view -h and then pipe this to htseq-count. Sorting the files prior to this conversion would cause errors.

The obvious first thing to check is that your gtf and genome match. Did you get them from the same place? Specifically, do the chromosome names in the gtf match the ones in the bam?

Good point..... I have just checked, and the gtf file contains no chromosomes actually. After running grep chr over the file it returns nothing. What is the reason for this? And how would I find one that works with the BAM files provided?

"grep chr" won't find anything if the chromosome names are bare numbers.  The question is, does your bam have bare numbers as chromosome names?

@swbarnes2 I converted your answer to a comment since it wasn't giving an answer but only requesting more information from the OP. Please only post answers if they are actually providing an answer. If you need more details, leave a comment instead.

Can you show us a few lines of your gtf file? I am guessing you have 1 instead of chr1.

My .bam file chromosomes are formatted as such: chr1, chr2,....

If your BAM file is using UCSC chromosome names (chr1, chr2, etc.) and your GTF file is using Ensembl chromosome names (1, 2, etc.) then htseq-count won't know how to match your alignments to any of the genes. The simplest solution is to not use htseq-count, but instead to use featureCounts, which will handle the chromosome name difference transparently (it's one of the few programs that does) and is much faster anyway. In the future, download your GTF and fasta files from the same source, don't mix and match them.

Thanks for the reply! My BAM file does indeed use this chromosome name system.Thing is though, my lab wants to continue using htseq for consistency. Would this problem be fixed by editing the .gtf file chromosome names to chr1, chr2 etc. ?

Yes, you can edit the GTF, though if you want real consistency you might see if you were using the exact same fasta and GTF files before.

@h3ab74 you can change the names by simply running sed 's/^/chr/' file.gtf > file.fixed.gtf`` and then use the file.fixed.gtffor your analysis. However, ***and this is important***: make sure that the bam files and the gtf are both referring to the same genome as Devon said, and also make sure there are no blank links in your gtf file. The command will only work correctly if _all_ the lines in the file should start withchr`.

I ran a similar command to that one: sed /^[[:digit:]]/s/^/chr/ file.gtf > newfile.gtf. I checked out newfile.gtf and made sure the chromosome names were corrected; they're now chr1, chr2, etc. However, when I ran the original command with the new file I still got the same problem. Is it maybe an issue with my samtools view command?

Thank you all for your helpful comments, I have made it work. And @DevonRyan, good point, I made sure to use the same .gtf files!

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


text:
REMC SNP annotations for transcriptome prediction - epiXcan
I have been trying to use a pipeline for predicting gene expression for a target gene (Although you can do it for many) - the name of this method is epiXcan. Here is a link to the paper:
https://www.nature.com/articles/s41467-019-11874-7#Sec22
One of the very first steps in the pipeline is obtaining SNP 'annotations' from an epigenome reference consortium known as the REMC. In the methods paper they say:
"We first prepare eQTL statistics (computed with MatrixEQTL4) and SNP annotations (extracted from REMC https://egg2.wustl.edu/roadmap/web_portal/). For each eQTL tissue, we use the matched REMC tissue to
extract the corresponding annotations".
I have completed the first step of obtaining eQTL statistics. However obtaining SNP annotations is where I am lost.
In my case, I am only trying to predict the expression of one gene (DRD4) and so my region of interest is the pre-frontal cortex. If you visit the website they provide though, it doesn't seem as though annotating a large list of SNPs - at least in the way they are presenting - is possible.
Can anyone provide some insight or direction on how to proceed? How does one obtain SNP annotations from the REMC? I have tried emailing the authors of the paper and they have not been responsive.
Thank you,

After much research I found out how to obtain the annotations.
You must go to the REMC website: https://egg2.wustl.edu/roadmap/web_portal/
Go to the chromatin-state-learning tab, and choose the n-chromatin state model desired - the core has 15 states, then 18, then 25 (imputed states).
Following this, you must go to the download URL link and, you will find a lot of files, primarily with the following coding information: E0*.. This corresponds to a particular tissue/cell type and so you must choose the one relevant to your research question. The way they are correspond to each tissue/cell type can be found in the following paper:
https://www.nature.com/articles/s41467-019-11874-7#Sec22


text:
Convert 10X VDJ .bam file to .fastq
I asked a similar question in the past regarding the conversion of scRNASeq bam files to fastqs. I am familiar with the bamtofastq functionality in the cellranger pipeline. However, it seems to not work with bam files from my VDJ samples (i.e. T-cell and B-cell receptor sequences), returning the following error with no other output:

bamtofastq v1.4.1
WARNING: no @RG (read group) headers found in BAM file. Splitting data by the GEM well marked in the corrected barcode tag.
Reads without a corrected barcode will not appear in output FASTQs'

Now, I read this biostars post indicating that this could be due to the fact that the bam file was not generated by cellranger, but I find that very unlikely given that it worked fine with the equivalent GEX (gene expression) sample.
So my question: how do I convert this VDJ bam to fastq? Most of the documentation I've read works primarily with GEX data, rather than VDJ, and given that bamtofastq worked fine for the GEX, why not with the VDJ? Is there a way around this?
Thanks,
UPDATE: I was able to edit the header of the file to add an @RG tag, however this did not fix the issue, and I got the following error instead:
Unrecognized 10x BAM file. For BAM files produced by older pipelines, use one of the following flags:
--gemcode   BAM files created with GemCode data using Longranger 1.0 - 1.3
--lr20      BAM files created with Longranger 2.0 using Chromium Genome data
--cr11      BAM files created with Cell Ranger 1.0-1.1 using Single Cell 3' v1 data
I'm beginning to question if this is even possible with VDJ data ..
P.S. This is the output of samtools view -h file_vdj.bam | head on the original bam file:
@HD     VN:1.6  SO:coordinate
@SQ     SN:clonotype1_consensus_1       LN:524
@SQ     SN:clonotype1_consensus_2       LN:570
@SQ     SN:clonotype2_consensus_1       LN:485
@SQ     SN:clonotype2_consensus_2       LN:514
@SQ     SN:clonotype2_consensus_3       LN:507
@SQ     SN:clonotype2_consensus_4       LN:521
@SQ     SN:clonotype3_consensus_1       LN:563
@SQ     SN:clonotype3_consensus_2       LN:507
@SQ     SN:clonotype4_consensus_1       LN:524
@SQ     SN:clonotype5_consensus_1       LN:563
@SQ     SN:clonotype6_consensus_1       LN:719
@SQ     SN:clonotype7_consensus_1       LN:501
@SQ     SN:clonotype7_consensus_2       LN:630
@SQ     SN:clonotype7_consensus_3       LN:511
@SQ     SN:clonotype8_consensus_1       LN:549
@SQ     SN:clonotype8_consensus_2       LN:557
@SQ     SN:clonotype8_consensus_3       LN:495
@SQ     SN:clonotype9_consensus_1       LN:524
@SQ     SN:clonotype9_consensus_2       LN:524
@SQ     SN:clonotype9_consensus_3       LN:587
@SQ     SN:clonotype10_consensus_1      LN:580
@SQ     SN:clonotype11_consensus_1      LN:570
@SQ     SN:clonotype11_consensus_2      LN:574
@SQ     SN:clonotype12_consensus_1      LN:563
@SQ     SN:clonotype12_consensus_2      LN:517
@SQ     SN:clonotype12_consensus_3      LN:568
@SQ     SN:clonotype13_consensus_1      LN:526
@SQ     SN:clonotype13_consensus_2      LN:563
@SQ     SN:clonotype14_consensus_1      LN:538
@SQ     SN:clonotype15_consensus_1      LN:524
@SQ     SN:clonotype15_consensus_2      LN:524
@SQ     SN:clonotype15_consensus_3      LN:571
@SQ     SN:clonotype16_consensus_1      LN:500
@SQ     SN:clonotype16_consensus_2      LN:562
@SQ     SN:clonotype17_consensus_1      LN:570
@SQ     SN:clonotype17_consensus_2      LN:753
@SQ     SN:clonotype17_consensus_3      LN:586
@SQ     SN:clonotype18_consensus_1      LN:524
@SQ     SN:clonotype18_consensus_2      LN:538
@SQ     SN:clonotype19_consensus_1      LN:556
@SQ     SN:clonotype19_consensus_2      LN:560
@SQ     SN:clonotype20_consensus_1      LN:456
@SQ     SN:clonotype20_consensus_2      LN:568
@SQ     SN:clonotype21_consensus_1      LN:524
@SQ     SN:clonotype21_consensus_2      LN:538
@SQ     SN:clonotype22_consensus_1      LN:432
@SQ     SN:clonotype22_consensus_2      LN:448
@SQ     SN:clonotype22_consensus_3      LN:580
@SQ     SN:clonotype23_consensus_1      LN:539
@SQ     SN:clonotype23_consensus_2      LN:621
@SQ     SN:clonotype24_consensus_1      LN:552
@SQ     SN:clonotype24_consensus_2      LN:489
@SQ     SN:clonotype25_consensus_1      LN:502
@SQ     SN:clonotype25_consensus_2      LN:512
@SQ     SN:clonotype25_consensus_3      LN:469
@SQ     SN:clonotype26_consensus_1      LN:654
@SQ     SN:clonotype26_consensus_2      LN:519
@SQ     SN:clonotype27_consensus_1      LN:527
@SQ     SN:clonotype27_consensus_2      LN:504
@SQ     SN:clonotype28_consensus_1      LN:508
@SQ     SN:clonotype28_consensus_2      LN:589
@SQ     SN:clonotype29_consensus_1      LN:501
@SQ     SN:clonotype29_consensus_2      LN:523
@SQ     SN:clonotype29_consensus_3      LN:618
@SQ     SN:clonotype30_consensus_1      LN:495
@SQ     SN:clonotype30_consensus_2      LN:517
@SQ     SN:clonotype30_consensus_3      LN:546
@SQ     SN:clonotype31_consensus_1      LN:563
@SQ     SN:clonotype31_consensus_2      LN:563
@SQ     SN:clonotype31_consensus_3      LN:507
@SQ     SN:clonotype32_consensus_1      LN:486
@SQ     SN:clonotype32_consensus_2      LN:577
@SQ     SN:clonotype33_consensus_1      LN:521
@SQ     SN:clonotype33_consensus_2      LN:538
@SQ     SN:clonotype34_consensus_1      LN:552
@SQ     SN:clonotype34_consensus_2      LN:647
@SQ     SN:clonotype35_consensus_1      LN:501
@SQ     SN:clonotype35_consensus_2      LN:514
@SQ     SN:clonotype36_consensus_1      LN:499
@SQ     SN:clonotype36_consensus_2      LN:618
@SQ     SN:clonotype37_consensus_1      LN:656
@SQ     SN:clonotype37_consensus_2      LN:525
@SQ     SN:clonotype38_consensus_1      LN:524
@SQ     SN:clonotype38_consensus_2      LN:526
@SQ     SN:clonotype39_consensus_1      LN:705
@SQ     SN:clonotype39_consensus_2      LN:482
@SQ     SN:clonotype40_consensus_1      LN:431
@SQ     SN:clonotype40_consensus_2      LN:580
@SQ     SN:clonotype41_consensus_1      LN:481
@SQ     SN:clonotype41_consensus_2      LN:535
@SQ     SN:clonotype42_consensus_1      LN:509
@SQ     SN:clonotype42_consensus_2      LN:489
@SQ     SN:clonotype43_consensus_1      LN:663
@SQ     SN:clonotype43_consensus_2      LN:589
@SQ     SN:clonotype44_consensus_1      LN:540
@SQ     SN:clonotype44_consensus_2      LN:508
@SQ     SN:clonotype44_consensus_3      LN:508
@SQ     SN:clonotype45_consensus_1      LN:479
@SQ     SN:clonotype45_consensus_2      LN:501
@SQ     SN:clonotype45_consensus_3      LN:456
@SQ     SN:clonotype45_consensus_4      LN:476
@SQ     SN:clonotype46_consensus_1      LN:662
@SQ     SN:clonotype46_consensus_2      LN:506
@SQ     SN:clonotype47_consensus_1      LN:458
@SQ     SN:clonotype47_consensus_2      LN:529
@SQ     SN:clonotype48_consensus_1      LN:484
@SQ     SN:clonotype48_consensus_2      LN:547
@SQ     SN:clonotype49_consensus_1      LN:532
@SQ     SN:clonotype49_consensus_2      LN:526
@SQ     SN:clonotype49_consensus_3      LN:508
@SQ     SN:clonotype50_consensus_1      LN:481
@SQ     SN:clonotype50_consensus_2      LN:495
@SQ     SN:clonotype51_consensus_1      LN:502
@SQ     SN:clonotype52_consensus_1      LN:502
@SQ     SN:clonotype53_consensus_1      LN:502
@SQ     SN:clonotype54_consensus_1      LN:552
@SQ     SN:clonotype55_consensus_1      LN:552
@SQ     SN:clonotype55_consensus_2      LN:500
@SQ     SN:clonotype55_consensus_3      LN:489
@SQ     SN:clonotype55_consensus_4      LN:551
@SQ     SN:clonotype56_consensus_1      LN:489
@SQ     SN:clonotype57_consensus_1      LN:540
@SQ     SN:clonotype57_consensus_2      LN:542
@SQ     SN:clonotype58_consensus_1      LN:541
@SQ     SN:clonotype59_consensus_1      LN:541
@SQ     SN:clonotype60_consensus_1      LN:541
@SQ     SN:clonotype61_consensus_1      LN:541
@SQ     SN:clonotype62_consensus_1      LN:512
@SQ     SN:clonotype63_consensus_1      LN:440
@SQ     SN:clonotype64_consensus_1      LN:500
@SQ     SN:clonotype65_consensus_1      LN:505
@SQ     SN:clonotype65_consensus_2      LN:533
@SQ     SN:clonotype66_consensus_1      LN:505
@SQ     SN:clonotype66_consensus_2      LN:458
@SQ     SN:clonotype67_consensus_1      LN:495
@SQ     SN:clonotype68_consensus_1      LN:495
@SQ     SN:clonotype69_consensus_1      LN:657
@SQ     SN:clonotype70_consensus_1      LN:645
@SQ     SN:clonotype71_consensus_1      LN:645
@SQ     SN:clonotype72_consensus_1      LN:495
@SQ     SN:clonotype73_consensus_1      LN:492
@SQ     SN:clonotype74_consensus_1      LN:492
@SQ     SN:clonotype75_consensus_1      LN:535
@SQ     SN:clonotype76_consensus_1      LN:534
@SQ     SN:clonotype77_consensus_1      LN:471
@SQ     SN:clonotype78_consensus_1      LN:540
@SQ     SN:clonotype79_consensus_1      LN:540
@SQ     SN:clonotype80_consensus_1      LN:473
@SQ     SN:clonotype80_consensus_2      LN:458
@SQ     SN:clonotype80_consensus_3      LN:502
@SQ     SN:clonotype80_consensus_4      LN:515
@SQ     SN:clonotype81_consensus_1      LN:522
@SQ     SN:clonotype81_consensus_2      LN:560
@SQ     SN:clonotype82_consensus_1      LN:492
@SQ     SN:clonotype82_consensus_2      LN:436
@SQ     SN:clonotype83_consensus_1      LN:657
@SQ     SN:clonotype83_consensus_2      LN:544
@SQ     SN:clonotype84_consensus_1      LN:547
@SQ     SN:clonotype85_consensus_1      LN:547
@SQ     SN:clonotype86_consensus_1      LN:557
@SQ     SN:clonotype87_consensus_1      LN:562
@SQ     SN:clonotype88_consensus_1      LN:562
@SQ     SN:clonotype89_consensus_1      LN:562
@SQ     SN:clonotype90_consensus_1      LN:465
@SQ     SN:clonotype91_consensus_1      LN:504
@SQ     SN:clonotype92_consensus_1      LN:498
@SQ     SN:clonotype93_consensus_1      LN:716
@SQ     SN:clonotype93_consensus_2      LN:511
@SQ     SN:clonotype93_consensus_3      LN:471
@SQ     SN:clonotype94_consensus_1      LN:485
@SQ     SN:clonotype94_consensus_2      LN:529
@SQ     SN:clonotype95_consensus_1      LN:519
@SQ     SN:clonotype96_consensus_1      LN:507
@SQ     SN:clonotype97_consensus_1      LN:521
@SQ     SN:clonotype97_consensus_2      LN:529
@SQ     SN:clonotype98_consensus_1      LN:536
@SQ     SN:clonotype98_consensus_2      LN:577
@SQ     SN:clonotype99_consensus_1      LN:502
@SQ     SN:clonotype99_consensus_2      LN:502
@SQ     SN:clonotype100_consensus_1     LN:471
@SQ     SN:clonotype101_consensus_1     LN:580
@SQ     SN:clonotype102_consensus_1     LN:580
@SQ     SN:clonotype103_consensus_1     LN:580
@SQ     SN:clonotype104_consensus_1     LN:580
@SQ     SN:clonotype105_consensus_1     LN:553
@SQ     SN:clonotype106_consensus_1     LN:553
@SQ     SN:clonotype107_consensus_1     LN:490
@SQ     SN:clonotype108_consensus_1     LN:491
@SQ     SN:clonotype109_consensus_1     LN:517
@SQ     SN:clonotype110_consensus_1     LN:501
@SQ     SN:clonotype110_consensus_2     LN:508
@SQ     SN:clonotype111_consensus_1     LN:510
@SQ     SN:clonotype111_consensus_2     LN:453
@SQ     SN:clonotype112_consensus_1     LN:494
@SQ     SN:clonotype112_consensus_2     LN:473
@SQ     SN:clonotype113_consensus_1     LN:545
@SQ     SN:clonotype113_consensus_2     LN:508
@SQ     SN:clonotype114_consensus_1     LN:505
@SQ     SN:clonotype114_consensus_2     LN:514
@SQ     SN:clonotype115_consensus_1     LN:659
@SQ     SN:clonotype115_consensus_2     LN:681
@SQ     SN:clonotype116_consensus_1     LN:491
@SQ     SN:clonotype116_consensus_2     LN:621
@SQ     SN:clonotype117_consensus_1     LN:517
@SQ     SN:clonotype117_consensus_2     LN:573
@SQ     SN:clonotype118_consensus_1     LN:513
@SQ     SN:clonotype119_consensus_1     LN:535
@SQ     SN:clonotype119_consensus_2     LN:537
@SQ     SN:clonotype120_consensus_1     LN:497
@SQ     SN:clonotype120_consensus_2     LN:579
@SQ     SN:clonotype121_consensus_1     LN:507
@SQ     SN:clonotype121_consensus_2     LN:519
@SQ     SN:clonotype122_consensus_1     LN:548
@SQ     SN:clonotype122_consensus_2     LN:526
@SQ     SN:clonotype123_consensus_1     LN:572
@SQ     SN:clonotype123_consensus_2     LN:565
@SQ     SN:clonotype124_consensus_1     LN:449
@SQ     SN:clonotype125_consensus_1     LN:490
@SQ     SN:clonotype126_consensus_1     LN:494
@SQ     SN:clonotype127_consensus_1     LN:578
@SQ     SN:clonotype127_consensus_2     LN:574
@SQ     SN:clonotype128_consensus_1     LN:571
@SQ     SN:clonotype129_consensus_1     LN:526
@SQ     SN:clonotype129_consensus_2     LN:571
@PG     ID:cellranger
@PG     ID:samtools     PN:samtools     PP:cellranger   VN:1.12 CL:samtools sort -o /lustre03/project/6033481/4nanuq2/platform/novaseqRuns/230605_A00977_0647_AH53N7DSX7_2155NS01A/10x_cellcount_1/CD8VDJ_MPS12349777_F10_2/SC_VDJ_ASSEMBLER_CS/SC_MULTI_CORE/VDJ_T_CLONOTYPE_ASSIGNER/CLONOTYPE_ASSIGNER/WRITE_CONSENSUS_BAM/fork0/chnk0-u059a808644/files/consensus.sam.sorted /lustre03/project/6033481/4nanuq2/platform/novaseqRuns/230605_A00977_0647_AH53N7DSX7_2155NS01A/10x_cellcount_1/CD8VDJ_MPS12349777_F10_2/SC_VDJ_ASSEMBLER_CS/SC_MULTI_CORE/VDJ_T_CLONOTYPE_ASSIGNER/CLONOTYPE_ASSIGNER/WRITE_CONSENSUS_BAM/fork0/chnk0-u059a808644/files/consensus.sam
@PG     ID:samtools.1   PN:samtools     PP:samtools     VN:1.12 CL:samtools view -bS -o /lustre03/project/6033481/4nanuq2/platform/novaseqRuns/230605_A00977_0647_AH53N7DSX7_2155NS01A/10x_cellcount_1/CD8VDJ_MPS12349777_F10_2/SC_VDJ_ASSEMBLER_CS/SC_MULTI_CORE/VDJ_T_CLONOTYPE_ASSIGNER/CLONOTYPE_ASSIGNER/WRITE_CONSENSUS_BAM/fork0/chnk0-u059a808644/files/consensus.bam /lustre03/project/6033481/4nanuq2/platform/novaseqRuns/230605_A00977_0647_AH53N7DSX7_2155NS01A/10x_cellcount_1/CD8VDJ_MPS12349777_F10_2/SC_VDJ_ASSEMBLER_CS/SC_MULTI_CORE/VDJ_T_CLONOTYPE_ASSIGNER/CLONOTYPE_ASSIGNER/WRITE_CONSENSUS_BAM/fork0/chnk0-u059a808644/files/consensus.sam.sorted
@PG     ID:samtools.2   PN:samtools     PP:samtools.1   VN:1.20 CL:samtools view -H NS.2155.001.SI-TT-A1_i7---SI-TT-A1_i5.CD8_VDJ.bam

Sorry, needs a few more lines, might be a lot more names depending on how many chromosomes there are. samtools view -h file_vdj.bam | head -n <linecount>. I'm interested in either a read group tag, or the first aligned sequence.

Thanks for your responses. I just added a little more! @gringer

This header doesn't look complete. You had SN tags previously, but they're gone now. Can you please show the full header? I expect it'll be relevant because the header is mentioned in the warning message.

I added the full header now @gringer

After contacting 10X genomics directly, I received the following response:
Please note that this is not 10x supported and is an untested work around.
Basically, what you can do is insert information from the gene expression BAM header into the VDJ BAM header and 'trick' bamtofastq to take the VDJ BAM. Let me outline the steps I took using the public data from https://www.10xgenomics.com/resources/datasets/integrated-gex-and-vdj-analysis-of-connect-generated-library-from-mouse-pbm-cs-2-standard-6-0-1.
Transfer VDJ BAM header to textfile
samtools view -H SC5v2_mousePBMCs_1Kcells_Connect_single_channel_Multiplex_vdj_b_all_contig.bam > vdj_b_all_contig.header
Transfer @RG and @CO GEX BAM header lines to textfile
samtools view -H SC5v2_mousePBMCs_1Kcells_Connect_single_channel_SC5v2_mousePBMCs_1Kcells_Connect_single_channel_count_sample_alignments.bam | grep -v '@SQ' | grep -v '@HD' | grep -v '@PG' > gex-rg-co.header
Concatenate the two text files
cat vdj_b_all_contig.header gex-rg-co.header > fixed.header
Replace the VDJ BAM header with the fixed header
samtools reheader fixed.header SC5v2_mousePBMCs_1Kcells_Connect_single_channel_Multiplex_vdj_b_all_contig.bam > fixed.bam
View the @RG line id field
grep '@RG' gex-rg-co.header
This will show the @RG line as follows for this public data I am using.
We need the highlighted content in the ID: field in the next step.
Add RG tags to each read
samtools addreplacerg -R SC5v2_mousePBMCs_1Kcells_Connect_single_channel:0:1:H3T2JDSX2:4 -o readgroup.bam fixed.bam
Run 10x bamtofastq on the readgroup.bam /Users/Documents/tools/bamtofastq_macos readgroup.bam testbamtofastqonvdj
This gives us matched R1 and R2 reads in the testbamtofastqonvdj folder.
One final note, in the final FASTQ set, be sure the R1 and R2 reads are written in the same order, i.e. the readnames match as you go down both. If you will map the reads with an aligner that uses insert size metrics, then it is important that the reads are not sorted by any genomic coordinate but rather are randomly sorted (samtools collate) or queryname sorted (samtools sort -n).


text:
Gene Ranking - signal to noise ratio used in GSEA-P algorithm?
I'm looking at Broad Institute's orignal GSEA-P algorithm R script which I downloaded here: http://software.broadinstitute.org/gsea/downloads.jsp.
I'm trying to adapt their GSEA.1.0.R script to process datasets that have 1 gene expression profile as for class 1, and N gene expression profiles for class 2. The bottleneck seems to be the GSEA.MetaboliteRanking() method. I looked and it appears they are using "signal to noise ratio" (SNR) metric to rank genes in patient profiles.
I was wondering how the signal to noise is calculated from gene expression profiles?
Note: I am using/modifying their script because it's the only implementation of the GSEA algorithm that I could find where they expose their code.
i.e., what is the mathematical reasoning behind the following lines of code? Note, there are little to no comments in the code, but the following variables, from what I can tell, are:
A: Matrix of gene expression values (rows are genes, columns are samples)
class.labels: Phenotype of class distinction of interest. A vector of
              binary labels having first the 1s and then the 0s

P = class.labels1 * subset.mask
n1 = sum(P[,1])
M1 = A %*% P
M1 = M1/n1
A2 = A*A
S1 = A2 %*% P
S1 = S1/n1 - M1*M1
S1 = sqrt(abs((n1/(n1-1)) * S1))

P = class.labels2 * subset.mask
n2 = sum(P[,1])
M2 = A %*% P
M2 = M2/n2
A2 = A*A
S2 = A2 %*% P
S2 = S2/n2 - M2*M2
S2 = sqrt(abs((n2/(n2-1)) * S2))

The noise is comming from the subset.mask, which is created above in a loop with the number of permutations.
for (r in 1:nperm) { #L90
....
        subset.mask[, r] <- as.numeric(c(subset.class1, subset.class2)) # L107

So by multiplying the random selection of subsets by the expression we get the "noise". Later on line 251 we get the ratio of noise/signal by dividing the signal from the contrast and the signal by the noise.
Those transformation on the middle seem to come from normalizing the signal in order to be able to compare the resampling, but I don't fully understand why they do this. Also note that the R scripts are outdated and no longer maintained, so they can contain errors.


text:
Recommendation for a binary dataset for Computational Biology open access (as in a dataset with 1 and 0 that can be used to apply ML techniques?
I have a couple of modified ML methods that I want to test with biological binary data - specifically binary gene expression data - e.g. https://eprints.soton.ac.uk/69359/1/Thesis.pdf that might include a variable that could be used to be the target variable in ML classification. What I mean by binary data is the features are all either 0 or 1, so I expect a large dataset e.g this thesis. I am looking for any suggestions.

Please [edit] your question and make it a bit more specific What kind of "biological data" do you need? How large a dataset? What do you mean by "binary"? Are you referring to the format of the data or do you want something that can only have two states?

Kaggle has a cats v dogs dataset if you want to try your hand at biological image classification: https://www.kaggle.com/c/dogs-vs-cats
The DREAM challenge competitions are usually specifically oriented to solving biological classification or learning problems with ML: http://dreamchallenges.org/


text:
Error in PSI-BLAST when running psipred locally
I am very new to bioinformatics and Linux and am trying to run psipred locally on Ubuntu 18.04 but cannot get PSI-BLAST to work. I have downloaded BLAST legacy version 2.2.26 and unpacked it, uniref90 database and psipred version 4.01
I installed psipred according to instructions:
cd psipred/src
make
make install

And this seems to work fine
I used the following code to create a database, which proceeded without errors:
makeblastdb -in uniref90.fasta -dbtype prot -out uniref90

In runpsipred I changed the directories:
# The name of the BLAST data bank
set dbname = /home/hilda/uniref90

# Where the NCBI programs have been installed
set ncbidir = /home/hilda/blast/bin

# Where the PSIPRED V4 programs have been installed
set execdir = ./bin

# Where the PSIPRED V4 data files have been installed
set datadir = ./data

I then try to run psipred:
tsch ./runpsipred example/example.fasta

And get the error
Running PSI-BLAST with sequence example/example.fasta ...
FATAL: Error whilst running blastpgp - script terminated!

I am not sure how to solve this as I do not understand what is going wrong. As I said I am new to Linux. I use the command tsch to run psipred (otherwise I get the error "command not found"), but not to install the database, could this be a problem? Or have I missed something in the setup process? Running psipred does create a copy of the example.fasta in the psipred folder named psitmp4286007f010.fasta, but that is all that seems to happen.
I would greatly appreciate any suggestions on how to make it work!
Edit:
I tried to run with BLAST+ (which seemed to work without tsch) and get the following error:
BLAST Database error: No alias or index file found for protein database [/home/hilda/uniref90] in search path [/home/hilda/psipred::]

I have also tried setting dbname to home/hilda/uniref90/uniref90 and got the same results. there is a uniref90.pal file in the uniref90 directory, which I understand is the alias file? Other forum answers seem to use pfilt to index (?) the database, but there is no pfilt in my psipred directory. Any ideas how to fix the problem?

Welcome to the site Hilda. I think that the supported version is now blast+ or something alike. Is there any special reason why you need the legacy 2.2.26? It will get easier to get support and trouble shot the problems if you install the supported version. In addition, it could be that the legacy version does not work with Ubuntu 18.04...

I use tcsh, but I would make a new attempt using bash.

I picked the BLAST version based on what seemed to be the most commonly used in forums and similar questions, and that I got the impression that psipred preferes legacy to +. Will definitely try downloading and running BLAST+ instead. Thank you!

I'm having the exact same problem. I'm using BLAST+ though and runpsipredplus script. Did you manage to make it run and how did you do it?

This error could be due to the fact that you are using legacy blast, have you tried using BLAST+ instead?
If you wish to use PSIPRED with BLAST+, then use the runpsipredplus script
here, rather than the normal runpsipred script.


text:
RNAseq data analysis
I am currently doing a Differential Gene Expression of RNAseq data in Lung Adenocarcinoma using TCGAbiolinks. In the data Preprocessing step TCGAanalyze_Normalisation, I am confused as to which method to select for normalization (geneLength or gcContent). What are the possible implications of both methods and how to decide which is the best suited model for the analysis?
Is there a better way to perform the DGE ?

The TCGAbiolinks reference manual provides a clear explanation of the two methods:

Normalization for RNA-Seq Numerical and graphical summaries of RNA-Seq
  read data. Withinlane normalization procedures to adjust for
  GC-content effect (or other gene-level effects) on read counts: loess
  robust local regression, global-scaling, and full-quantile
  normalization (Risso et al., 2011). Between-lane normalization
  procedures to adjust for distributional differences between lanes
  (e.g., sequencing depth): global-scaling and full-quantile
  normalization (Bullard et al., 2010).

Ideally, you would know if your samples spanned different lanes or not, which would indicate which you should use - "geneLength" if so, "gcContent" if not. If you're using TCGA data, that metadata is likely available somewhere.


text:
Retrieve RNA sequencing data for human p53 colon cancer cell lines
I want to download the next generation RNA sequencing data in different cell lines. I want the data (alongwith normal) for p53 wild type, p53 knock-out and p53 null cell lines (human colon cancer). From where do I download this data? In fastq format if possible.
Currently, I am going through SRA, to look for the datasets, using their search engine. But the process is extremely tedious.
Also can someone please hint about how to download the data from SRA, in case that's the only way to get data. I'm getting very confused.

Welcome to the site. Could you help us explaining if you could manage to download any dataset? Also, what kind of sequencing data do you want: single-cell RNa-seq, long read or bulk RNA-seq? Also if you could clarify if you have some preference for a programming language to automatize the process it would help too

Hello, What exactly you find tidious? If you search for "p53 AND "Homo sapiens"[orgn:__txid9606]" you get plenty of hits for NULL and WT. For example reads liked to this project. If you want fastq files, you can pull them from ftp server once you know IDs, here is an example: ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR835/SRR835937/

Hello, thank you all for such an immediate response. Ilrs, I want RNA seq data obtained from p53 WT, NULL and KO colon cancer cell lines, without any other modifications in the cell. I'm not bothered about programming language for further analysis now. I am actually new to this, and I don't know if I am collecting the right kind of data or not.  Kamil, these are shPFKFB4 silenced. I want normal vs experiment data from same experiment if possible. I might be making very silly mistakes, I apologise for the same. Thanks for your help.

Downloading from ENA is simpler: https://www.ebi.ac.uk/ena

I suggest you to try difference resources, including ArrayExpress.
A simple search at ArrayExpress returns multiple experiments which might be of your interest. Just a pointer, ArrayExpress is also synced with GEO.

I agree with Kamil on the downloading method.
Moreover, SRA will give you your file in the raw format which is usually a Binary Alignment Map (BAM) or Sequence Alignment Map (SAM)
Download SRA toolkit available online, to convert your raw reads to fastq format. Follow the steps using command line. It is a very simple command to convert your SRA file to fastq which uses the fastq-dump command.
Further reading: https://www.ncbi.nlm.nih.gov/books/NBK158900/
Hope it helps!

Hello Pawan, my main concern here it to pinpoint the data I require from the database. I am not able to find the data I require. I'm aware about fastq-dump. Thank you for your efforts.

I understand! Since you had asked if it was possible to download in fastq format, i had answered accordingly. Cheers

Thanks. Actually, when I'm using fastq-dump, all I see upon execution of command is a bunch of .cache files.Also, as huge files in GB scale take time to download, I would like to know if there is a way where we can track the progress of download?

Can you specify the size of your file?

Last file I tried downloading  was roughly 12GB. Max download speed on my network is 550kbps.

Well, downloading at this speed will take over 6 hours.


text:
set -euo pipefail in Snakemake wrapper
I have snakemake script doing fastqc with a wrapper but it errors out saying
Command 'set -euo pipefail;  /home/mi/ppp/miniconda3/envs/meta/bin/python3.9 /workflow/.snakemake/scripts/tmpsl695gpr.wrapper.py' returned non-zero exit status 1.
  File "/workflow/rules/fastqc.smk", line 32, in __rule_fastq1
  File "/home/mi/ppp/miniconda3/envs/meta/lib/python3.9/concurrent/futures/thread.py", line 58, in run

my rule looks like this
rule fastq1:
    input:
        "../resources/samples/{sample}_1.fq.gz"
    output:
        html = "../results/html/{sample}_1.fq.html",
        zip = "../results/zip/{sample}_1.fq.zip"
    conda:
        "../envs/meta.yaml"
    log:
        "../results/logs/fastq_raw/fastq_raw_{sample}.log"
    threads:
        2
    wrapper:
        "v1.5.0/bio/fastqc"

I have a feeling this used to work before and my partner says it works fine on his machine but I can't find the issue.
I'm also doing this on a clean conda/mamba install and env.
EDIT:
Snakemake version is 7.7.0
I "solved" this by installing fastqc with conda... I mean am I stupid? Is it a requirement for a wrapper to also have the package installed? I don't think this was mentioned on the Snakemake website
(Im not sure wht mr communuty bot is complaining about the problem is the error message and me not being able to run the code..)
Thanks for reading.

Please include Snakemake version in question body.

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.

Note that set -euo pipefail is something Snakemake adds by default to your shell commands. The fact that it appears in the error message is usually not relevant. The error likely has to do with something else in the command that the wrapper calls.


text:
Demultiplex nanopore reads with custom barcodes
We have a problem trying to demultiplex MinION sequences with custom barcodes. Do you have any software recommendations we can try for demultiplexing or how to demultiplex these custom barcodes with Albacore? We have tried using albacore but it only recognizes barcodes provided by Nanopore Technologies.

I've written some Perl code for demultiplexing nanopore reads using LAST. A workflow demonstrating the use of this script (and how it works under the hood) can be found here:
https://dx.doi.org/10.17504/protocols.io.14egnxw4zl5d
You feed the script a barcode file and input fastq file, and it creates a trained mismatch matrix, aligns to the barcodes, then demultiplexes the reads into separate files for each barcode:
./fastq-dental.pl -barcode barcode_full_PBK004.fa -mat bc.mat reads_all.fastq.gz

Example output:
         9 BC01
        39 BC02
         1 BC02_BC07
        16 BC03
        38 BC04
         2 BC04_BC06
         1 BC04_BC07
         1 BC04_BC09
       643 BC05
         2 BC05_BC07
    493175 BC06
       115 BC06_BC06
         2 BC06_BC07
    663161 BC07
         1 BC07_BC05
       117 BC07_BC07
         1 BC07_BC09
       106 BC08
         1 BC08_BC06
        75 BC09
        72 BC10
         1 BC10_BC07
        41 BC11
        44 BC12
         1 BC12_BC06
     13535 BCchim
    237796 BCnoadapt
       152 BCnone
        90 RB12A

While the approach designed around using the ONT barcodes, the barcode fasta file can be replaced by another file [using barcode ID in the header] without change the protocol.
FWIW, there is also now some support for detecting strand direction (e.g. using the strand switch primer provided by ONT), but I haven't yet got it working properly doing both demultiplexing and strand orientation at the same time.

Thank you! I am using minibar, works great and very easy to add custom barcodes.
L.Y.

There are two other ways in which you could demultiplex reads if guppy doesn't have the config file for the kit you are using.

Minibar: This works for dual index barcodes specifically
Porechop: This is intended for the ONT sets, but you can also replace them with your own. However, it is designed to work with single index barcodes. That means, you will have to run it once with your forward indices, then take that output and run again with the reverse, and then do some gymnastics to combine it all together

What do you mean with run it twice, with forward and reverse indices?


text:
Repetition-free Longest Common Subsequence as similarity measure
I'm starting a thesis in Computer Science with topic the RFLCS (Repetition-free Longest Common Subsequence). I know very little about Bioinformatics problems and Genomics. I would like to understand the motivations behind this problem.
More precisely, I'd like to understand why one need a similarity measure between genomes, and how the RFLCS is an useful measure.
[Adi et al.] introduce the problem in this article.

What is the RFLCS? PLease [edit] your question and give us some more detail so we can understand exactly what you're asking. However, I am afraid that this will probably be too broad to answer (see the "too broad" section here: https://bioinformatics.stackexchange.com/help/closed-questions).

Done! I narrowed the field, hope enough.

Still a bit obscure. I think the question is there to sharpen your coding skills. Its fairly easy to solve. Although it could be solved as a 5% tail of a Poisson distribution (I doubt that is what its about).

I found the paper introducing it and I think you should mention it in the question. I think that your supervisor is the best person to talk to regarding motivations of the problem. If there is anything you would like to get help with from the community, you need to be more concrete. "Why do I do my thesis" is way too broad.

A brief literature review can give you the answers.
This preprint gives the definition of RFLCS more clearly:

A repetition free Longest Common Subsequence (LCS) of two sequences x and y is an LCS of x and y where each symbol may appear at most once.

Adi et al proposed to solve the genome rearrangement problem raised in this paper. If you want to know more, follow the references and citations of these papers and read them.
Basically, Adi et al think the length of RFLCS measures the distance between two genomes. The longer the more similar. I am not sure if this is a good measurement, though.

Now I see. If shifts in LCS were solely the consequence of rearrangements .. yes, but we all know thats not true. Transposons move, microsatellites, minisatellites and DNA virus integration (e.g. HPV) will all be confounders,

Thanks for the answer. I read several papers, but no one explains why RFLCS should be a good measure. Seems nobody have assessed its effectivness.

@DanieleCuomo I don't think RFLCS will be effective, tbh.


text:
GTF from consensus sequence
I am new to using bioinformatic tools and I was hoping this community could help clear something up for me. I need to generate a gtf file. My data are a set of complete HA genes for influenza B viruses in fasta format. From reading through forums and other internet searches, my understanding is that fasta cannot be transformed to gtf or gff because these last two provide location information which fasta does not (unless available in the header). I saw someone suggest mapping sequences to a reference sequence and generating the gtf from that mapping.
My question is, would it make sense to use the consensus sequence selected for alignment as this reference sequence? The alignment algorithm is a simple majority rule. Because all my data are complete genes, I would think the mapping of locations (start codon and stop codon) would be the same no matter which sequence in my data set I choose to use. Is this not so?
I am looking at viral evolution in the HA segment and my end goal is to calculate the dS/dN ratio for a set influenza B viruses collected during a particular season. The program I want to use requires a gtf, but likely other options would only require the reference genome. Would it make more sense to use the strain used as the vaccine candidate as the reference genome?
I appreciate any help you can provide me!

Welcome to the site. I am not used to work with gtf, but you'll probably need a reference genome in order to align your sequences to the genome. Do you have a reference genome for the influenza B virus? What I know is that viruses change quite a lot and probably there are several genomes available. You seem familiar with alignments but you way that it is a matter of "simple majority rule", some algorithms are designed for short sequences or faster while some other are optimized for more accurate mapping, so it is not just a majority rule, but others will answer you with more detail than I can.

I think the answer is no. Let me check with my colleague who has more expertise in gtf or gff, before I cite my reasons and suggested way forward. Are you working within or associated with Derek Smith's lab (Cambridge)?

In summary, yes to this "I saw someone suggest mapping sequences to a reference sequence and generating the gtf from that mapping." and no to this "My question is, would it make sense to use the consensus sequence selected for alignment as this reference sequence?"

@MichaelG. Thank you for looking into this for me! No, I am not working within or associated with Derek Smith's lab in any way. Could you expand on why using the consensus sequence does not make sense?

I would still like to understand why using the consensus sequence does not make sense, but my other way forward is to use the sequence for the strain used as the vaccine candidate for that particular season. Does this make more sense? Because influenza has such a high mutation rate, it would seem that having a "standard" reference genome would not apply as it does for human genomes.

Is making a gtf file really your end goal?  I have the feeling that you are trying to do something else, and think that you need a gtf, but maybe that's not the case.

@swbarnes2 My end goal is to calculate the dS/dN ratio and the program I want to use requires a gtf. I imagine other software options would not require a gtf and instead would ask for a reference genome.

Welcome @Lindsey! It looks like these comments are getting somewhere. Please use the [edit] button above to update your original question to include a bit more about what you're trying to accomplish. It looks like you've sequenced some flu strains and would like to look at the evolution of HA genes?

I was wanting more time, but okay here's what I'm thinking,

Would it make more sense to use the strain used as the vaccine
  candidate as the reference genome?

Neither, because the GTF format is based on an HA protein structure (likely a crystal structure but could be cryoEM). This is not the same as an alignment position within a surface antigen, because most alignments of pathogen surface antigen genes will contain indels, e.g. between the GTF and your influenza epidemic. If HA doesn't have any indels fine, but I would be surprised.
The indels will cause artifactual shifts in the GTF domain boundaries unless they are corrected. Put simply an alignment position (particularly for pathogen surface antigen proteins) is rarely the same as the amino acid position of GFT because a protein structure of a surface antigen doesn't comprise indels in its native form.
If you do this work in detail, i.e. shifting between protein structure and phylogeny, you always write scripts which convert between aligned positions and "reference" (crystal/ cryoEM structure) positions. Its a common way of thinking.
So in summary, your reference sequence is the GTF sequence because that is the basis of the structure<-> alignment conversion

For positive selection analysis it does not matter what the reference is, all that matters is the correct rooting. For demarcation of domains within HA the GFT sequence needs to be at the top of the alignment 1) because it helps easily identify whether it contains indels 2) structure programs will expect the reference sequence to be there. If you are using the vaccine in a positive selection analysis ... you need to think hard about the question you are asking and the coverage of the vaccine. Epidemiologically speaking thats complicated because the passage of the vaccine will be an important source of selection and you will mistake that as the wild-type "escaping" the vaccine.

What I was suggesting was using the sequence for the HA segment of the vaccine candidate as the reference sequence used for the mapping and subsequent generation of the gtf. I'm not sure I made that clear. As you say, the HAs of the epidemic are likely to have a lot of indels, but still, some sequence must be considered the reference on which to align. In viruses with high mutation rates, what would make a good reference sequence?

I've amended above

So if I am understanding you correctly, I could use either the consensus sequence or the vaccine candidate HA sequence (or any sequence) as the reference in a positive selection analysis as long as I have the correct rooting. I'm not sure if I fully understood the last part of what you amended, but my set of HA sequences come from influenza cases, none are from the vaccine.

Hi Lindsey, .... firstly NEVER use the consensus sequence in a positive selection analysis - even if we disagree about the epidemiological focus of your study - NEVER use consequence in a phylogenetics analysis. However, yes the maximum likelihood approach to positive selection (the only one worth bothering about), the correct root is important. In other positive selection calculations (distance), it is will not make much difference.

"Transforming" Fasta files to GTF or GFF3 files is a common request, as are similar tasks such as "converting" BAM files to VCF files. But these are underdefined questions, and "transform" or "convert" isn't really the right way to think about the task.
As far as your particular request is concerned, Fasta files and GTF files contain fundamentally different types data.

Fasta files contain sequences, which could be chromosomes, assembled contigs, polypeptides, transcripts, or other biological sequences.
GTF (and GFF3 and BED) files contain intervals, indicating that some feature of interest occurs at such-and-such position on such-and-such sequence. A hypothetical example: a transcription factor binding site occurs on chromosome 2 from position 40011 to 40017.

So you are correct: it's not really helpful to think of your task in terms of "transforming" Fasta to GTF. The process of creating a GTF/GFF3/BED file is often called annotation. There are many different methods to annotate many different types of genomic features: transposable elements, segmental duplications, protein-coding genes, non-coding RNAs, and so on.
As a result, it's meaningless to say "I need a GTF file" without clearly specifying what that GTF file will annotate. (It's kind of like saying "I need an XML file" or "I need an Excel spreadsheet" without indicating what the contents of those files would be.) The things you need to clarify are:

What will the intervals in my GTF file represent?
What precise sequences will serve as the reference and provide interval coordinates for the GTF file?
How are those intervals to be determined?

It sounds like 3. is a big part of your question. If I'm reading things correctly, you want to annotate a particular class of genes on an influenza genome. Is that correct? Did you perform de novo sequencing of influenza viruses, and now hope to characterize variation in their HA genes? Stating what you hope to get out of your analysis would be helpful in finding a good way forward.

Thank you for replying! I have updated my question to be more informative of what I want to accomplish.

If I was doing this analysis, I'd want to assess selection signal against domain boundaries. In a selection calculation its one approach to demarcate the loci. For example the hydrophobic "stem" (I don't know HA terminology) will not be of huge interest in this particular question, so could be delineated to provide a control

Since you know that your sequences in your fastas are your genes of interest, why not just dummy up a gtf where every sequence is a gene?
The point of a gtf is to tell software which genomic regions matter for a certain procedure, and which regions do not.  But if you already know that your sequences are pared down to genes, then you know that the whole thing is what's relevant.


text:
How to use analysis output in my pipeline
I am running an analysis and I want to edit the output that I received from that analysis. However I cannot define the output name in my pipeline, because after running the analysis it is asking me to name the output file. How can I include that output file in my pipeline ?

Welcome to Bioinformatics.StackExchange! I would recommend reverting the edit you just made. The question is impossible to answer generically, since many tools have different ways to specify output filenames. My answer is only applicable to pepstats, like your original question asked.

I assume you are referring to the pepstats program from EMBOSS. The program has a -outfile flag that allows you to specify your desired name for the output file when you're not running the program interactively.


text:
I need some tips and suggestions for further analysis of NGS expression data (log2cpm)
I am a PhD student who inherited some log2cpm data of expression data from bulk kidney tissue from a UUO(unilateral urethral obstruction) experiment that tests a new drug. The sample material consists of:

6 x Ligated Kidney (Untreated
6 x Ligated Kideny (Treated)
3 x Unligated Kidney (Untreated)
3 x Unligated Kidney (Treated)

The previous study showed that this drug was effective against fribrosis, so the aim of my study is to investigate how this drug affects inflammation and mitochondrial function in the diseased kidney.
So I have filtered my original dataset based on MitoCarta v2 DB, so investigate how the drug affected mitochondrial function. I afterwards run a PCA analysis on log2cpm data which revealed the following:

The following figure shows PC1(84.51%) and PC2(4.69%).
Red (Ligated, Untreated), purple (Ligated, Treated), green(healthy kidney, not treated), orange (Healthy kidney, treated).
I interpret this as there is a very small difference between treated and not treated ligated kidneys, however there is a small difference.
I have just started doing bioinformatics and I am very unsure which approach I should choose further, however I have though maybe do a highly variable genes analysis (HVG) see link, isolate the 500 most variable genes and run a heatmap with hieraichal clustering. I was also thinking of looking at the loadings in PC2, and take this further. The truth is that I am very unsure what is the right approach, and my research group does not have bioinformatician I can ask anymore.
I have the following questions:

What kind of analysis you guys would suggest I use to investigate
further (Hieraichal clustering, HVG etc?)
Is there any alternative preprocessing I could use before I run the PCA, to make the analysis more accurate?
Are there any packages that you can suggest?
Any other tips?

Any other tips?
I asked this same question on Biostars as well
https://www.biostars.org/p/367722/#367730

Welcome to the site. What is your expected result? A list of genes to then validate experimentally with an immuno* or some list of genes for prediction of status/response... What other information do you know about those samples: Could it be that some samples are before and after treatment? Have you considered analysing the differentially expressed genes on the whole dataset (without filtering first for mitochondrial related genes)?

It looks like you'd be best served with a standard differential expression analysis using a factorial design. You should get the raw counts data for that though.

Hey,
the  Animals were sacrificed on day 7 or day 15 after the surgery, and treatment started immediatley after surgery. Drug was given immediatley after surgery.
The reason I want to investigate mitochondrial and immune genes, is that I previously did a GSEA with TopGO package of the whole transcriptome based on logfc-values and q-values (down and upregulated genes isolated, qval-cutoff 0.01, (+/-)fc 1.5) of Ligated(untreated vs. healthy kidneys), and Ligated(treated) vs. Ligated(untreated) kidneys. see figure. It showed that mosltey mitochondrial and inflammatin was affected

Now I want to investigate a bit further what molecular pathways the drug affects and what constitute the improvement.

@NewbieCoder Did you include the time in the design matrix when analysing differential expressed genes? Also if you are interested in those pathways is different than being interested in the genes expressed in the mithocondria as it appeared first. Could you [edit] the question and clarify what is your question?

Hey, I can add some information tomorrow. I was not at work the past days :) Then I can clarify exactly what I need to do. If you have any tips, I would appreciate it greatly :)

Since you indicate already having performed differential expression and using gene set enrichment the next step is to look at what changes due to drug treatment in the pathways of interest while having a decent understanding of those pathways. Things like IPA or pathview (from Bioconductor) can aid in this, but at the end of the day you need to have a decent understanding of the pathways yourself to judge what seems to be like the best targets for future work.


text:
How to annotate gene length to a list of gene symbols using UCSC data?
I have a list of HGNC gene symbols, I am looking to get the gene length of each gene. Although I also describe these genes with lots of UCSC datasets as features, so I am wondering if there is a dataset in UCSC I can use to also get gene length from?
I've been looking in the data that is downloadable from UCSC table browser (I've been aiming to find start and ends for each gene to subtract to get gene length) but there are a lot of files and I'm not sure which dataset to take from which will also match to my HGNC gene symbols.

You could do a few command-line operations to answer this question. This assumes the use of hg38 assembly.
First, get a list of genes from GENCODE:
$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "gene"' - \
    | convert2bed -i gff - \
    > genes.bed

Then use grep to filter those genes with your list of HGNC symbols:
$ grep -wFf hgnc_symbols.txt genes.bed > filtered_genes.bed

You can run this through awk to get lengths:
$ awk -vFS="\t" '{ print $3-$2 }' filtered_genes.bed > filtered_gene_lengths.txt

Thank you for this, this looks great. I should've said I am using hg19, my own fault, but I can figure it out from here so thank you :)


text:
How to compact variant data to their genes without bias?
I have a dataset of genes I am trying to collect data on from public databases, to use as features in machine learning. I am trying to take some features from UCSC genome browser (e.g. number of CpG islands per gene, number of DNase clusters per gene, regulatory enrichment scores etc.) however I am not sure how to control for bias where a gene that is larger in length - and so will then have more CpG islands or higher regulatory enrichment scores simply due to gene length.
Is there a way to correct for gene length when taking/condensing variant data to individual genes?
Edit:
For reference, my machine learning model aims to predict whether a gene is the most likely to be causal for a disease (out of all the genes given to the model). The model will score the genes as a regression classification between 0 to 1 (0 being least likely to cause disease and 1 being most likely to cause disease). I plan to later further investigate the genes with the highest scores.
The model uses a variety of multi-omic features (e.g. GTEx gene expression the genes have for many tissues, GWAScatalog data, gene intolerance scores, protein-protein interaction data, drug interaction data, phenotypic scores etc.). However, I am missing epigenetic data to describe my genes so I've been looking to collect based on UCSC's variant data (CpG islands, histone modifications, DNase clusters) - however this leads to my gene length problem when I am trying to reliably take data from the variant level.
I've been plotting my features and gene length, and seen that the UCSC epigenetic data does correlate with having a larger gene length if there is a higher count of regulatory sites (0.8 r2 for some), and so this is what I'm looking to correct.

Does regulatory enrichment score correspond to the number of regulatory SNPs / other regulatory elements per gene?

I believe so, a score is given per variant in the public database. I am then finding all those scores per gene in my gene list and usually just picking the highest score to be the value I use for that gene, only problem being a longer gene will have more opportunity for variants and more scores. Although, I am more focused on just creating counts of different regulatory sites per gene (e.g. number of CpG islands) and using that count as a feature input for machine learning (so I can not use the scores at all) but I am still stuck with the gene length issue.

What do you want to do with your machine learning? (what's correct depends on what are your goals) Have you looked at the way cqn package at Bioconductor corrects for gene length?

I am looking to use machine learning to predict whether a gene is most likely to cause a disease (model scores the gene with a 0-1 likelihood, predicting that score based on multiomic features I collect to describe the gene). I have looked at cqn but its starting point is RNA-seq data which I don't have. I only have gene symbols and the variants/SNPs within those genes that have come from GWAS to associate to the disease. Thank you for your comments, I'll add this info into my queston for clarity.

Its very easy, just let the ML sort this out for you and that is its advantage, You're thinking of GLM style calculation where you pre-screen the data with bivariate plots, where there needs to be nice Q-Q plots and low residual.
For ML simply include the gene length as one of your parameters along with CpG etc ... and the ML regression analysis SVC, lasso, ridge, random forest will figure the relationship out between gene length and CpG. You do zero, the ML does everything, hence from a statistical point of view purists object because you don't know the relationshiop the ML has deduced between the variables, but you will get regression weights for non-DNN stuff, which will give you some idea of the impact of length.
There is the issue of transformations and that can be complicated, but I'd try untransformed data first. The only disadvantage of this approach is the user will have to input the gene size when they want to check out your training algorithm.

Thank you for this, this is ideal for me. I'd give you more than 100 points if I could (when it lets me I will award the 100 points). I have a final question for you if possible: when talking about this it was originally suggested I look into statistical tests (I have a biology background so have been unsure) like Fisher’s Exact Test or Chi2. I didn't find anything I thought would easily fit my data - is the GLM style calculation you first mention the only type of statistical testing I could do for this? I can't find any commonly used tests which don't start with RNA-seq data.

My advice is to stick with ML but do this well, and you don't need a 'fit' anything - you definitely don't need to worry about GLM (ML is more trendy in any case). Geting the ML working to a good standard can be tricky and I'll discuss the issues above shortly.

Thank you for this, I will stick with ML. Currently I use xgboost and view the feature importance, is there a way for me to prove xgboost is regulating how it handles certain features in relation to the gene length feature? I'm trying to finding examples of this but potentially this is a new question for me to ask/research.

I thought xgboost just handles gradient descent for the ML regression. Are you integrating this with sci-kit learn? Using different methods is important here and further discussion is a good idea.

I use xgboost integrated with scikit-learn, would it be better done with xgboost functions only? I also do an initial benchmarking against xgboost, random forest, GBM, decision tree, k-nearest neighbour and logistic regression (elastic net and lasso). Although xgboost is always the top performer (with 10-fold nested cross validation) when I've tried different things (like different feature selection methods or hyperparameter tuning methods).

No using sci-kit learn is essential.. What about SVM? What is your accuracy? It is unusual when a single approach easily outperforms everything else in my experience.

I've got an r2 of 0.8 with xgboost and 0.75-0.6 r2 with all the other treebased methods. k-nearest is around 0.6 and the logistic regressions are ~0.2 r2. I did a while ago try SVM but found it had extremely low r2 no matter what I did (but I plan to revisit this as I have new/bigger data and hyperparameter tuning). Also I look at other metrics as I know r2 isn't the best, like the mean square error (from 0.16 - 0.18 for the tree-based methods for example). Thanks again for your help as I know this is moving away from the question, regulating gene length has been one of my last puzzle pieces

I'll have another look at your comment later, but I think that r2 is not 'accuracy'. Accuracy is the number of correct positives over the total number of predictions performed on the 'test' split of the data. An accuracy of 0.8 for a biological/genetic system is good enough.and if you have a method to achieve this I think most would suggest this is your result. There may not be universal consensus on this though because some argue that the 'consensus' of all methods is the result that matters (unless it is a DNN).

Thank you I'll look into what I can do with getting an accuracy. My only concern is with the regression say the model predicts a gene is scored at 0.801 when the true training score was 0.8 exactly - do I count that as a correct positive or would the model need to be exactly aligned with the real score?

The 0.8 has the look and feel of an accurary score. I wouldn't worry about the 0.801 issue. Accuracy is the primary index used in ML


text:
Querying metadata (GDC) using a filter
I have been trying to access the GDC metadata with python through the filter below. However, this code I have been using always yields a blank tsv file. Any help on fixing the code is appreciated!
import requests
import json

cases_endpt = 'https://api.gdc.cancer.gov/cases'

# The 'fields' parameter is passed as a comma-separated string of single names.
fields = [
    "sample_id",
    "sample_uuid",
    "sample_type",
    "sample_type_id",
    "tissue_type",
    "tumor_code",
    "tumor_code_id",
    "oct_embedded",
    "shortest_dimension",
    "intermediate dimension",
    "longest dimension",
    "is_ffpe",
    "pathology_report_uuid",
    "tumor_descriptor",
    "current weight",
    "initial weight",
    "composition",
    "time_between_clamping_and_freezing",
    "time_between_excision_and_freezing",
    "days_to_sample_procurement",
    "freezing_method",
    "preservation_method",
    "days_to_collection",
    "portions"
    ]

    fields = ','.join(fields)

filters = {
    "op": "in",
    "content": {
        "field": "project_id",
        "value": ["TCGA-BRCA"]
    }
 }

params = {
    "filters": json.dumps(filters),
    "fields": fields,
    "format": "TSV",
    "size": "10000"
    }

response = requests.get(cases_endpt, params=params)

file = open("query10.tsv", "w")
file.write(response.text)
file.close()

I think you are using an outdated version of the API, you need to change the "field" with cases.project.project_id
import requests
import json

cases_endpt = "https://api.gdc.cancer.gov/cases"

fields = [
    "sample_id",
    "sample_uuid",
    "sample_type",
    "sample_type_id",
    "tissue_type",
    "tumor_code",
    "tumor_code_id",
    "oct_embedded",
    "shortest_dimension",
    "intermediate dimension",
    "longest dimension",
    "is_ffpe",
    "pathology_report_uuid",
    "tumor_descriptor",
    "current weight",
    "initial weight",
    "composition",
    "time_between_clamping_and_freezing",
    "time_between_excision_and_freezing",
    "days_to_sample_procurement",
    "freezing_method",
    "preservation_method",
    "days_to_collection",
    "portions",
]

fields = ",".join(fields)

filters = {
    "op": "in",
    "content": {"field": "cases.project.project_id",
                "value": ["TCGA-BRCA"]},
}

params = {"filters": json.dumps(filters),
          "format": "TSV",
          "size": "10000"}

response = requests.get(cases_endpt, params=params)
content = str(response.content, "utf-8")

with open("query10.tsv", "w") as file:
    file.write(content)

And if you use with open, you don't need to use close.


text:
Trimmomatic trims most of the reverse strand
I have a paired-end data. It is from small part of human genome not a WGS or WES. I use BWA for alignment and do variant calling. Since I had some false positives I wanted to do trimming.
I used trimmomatic with command:
trimmomatic PE -phred33 -threads 2 read1.fastq.gz read2.fastq.gz read1_trimmed.fastq read2_trimmed.fastq LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15`

I do not know why but I have this output:
Input Read Pairs: 1834513
Both Surviving: 1822034 (99.32%)
Forward Only Surviving: 5740 (0.31%)
Reverse Only Surviving: 417 (0.02%)
Dropped: 6322 (0.34%)
TrimmomaticPE: Completed successfully

I do not have any adapters since the are removed in bcl2fastq step. In FASTQ quality control I have seen that quality of reverse and forward reads are similar in quality but trimmomatic trimmed many reads in reverse file. I cannot do alignment with these files, BWA gives a warning that says:
[bseq_read] the 2nd file has fewer sequences.
[process] read 155012 sequences (16693324 bp)...
[bseq_read] the 2nd file has fewer sequences.

Is there anyone who faced a problem like this before?
Thank you in advance

Maybe you need to remove the unpaired reads:
trimmomatic PE -phred33 -threads 2 \
    read1.fastq.gz read2.fastq.gz \
    read1_trimmed.fastq.gz read1_unique.fastq.gz \
    read2_trimmed.fastq.gz read2_unique.fastq.gz \
    LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15

Then you can try with bwa.