text:
Boxplot of read length distributions from multiple samples
Newbie R user here. I want to create a plot with boxplots of read length distribution on the Y-axis, and samples on the X-axis -- like this image from doi: 10.1038/srep16498
However, I'm not sure how to plot my data in ggplot2 to make this. My data is currently in this format:
Length Sample1 Sample2
20 10 5
30 15 223
40 16 4923
50 293 239
60 3290 23
70 248 6
80 23 1
90 21 0
100 5 0

Length is the read length bin, with the sample columns containing the number of reads of that length. I've tried using ggplot2, but don't think I can do it with the data in its current format.
Thanks!
Edit: Here are the number of reads for each individual read length (converted to long format):
"Readlength" "name" "value"
"1" 25 "Sample1" 457
"2" 25 "Sample2" 580
"3" 26 "Sample1" 603
"4" 26 "Sample2" 648
"5" 27 "Sample1" 599
"6" 27 "Sample2" 715
"7" 28 "Sample1" 686
"8" 28 "Sample2" 914
"9" 29 "Sample1" 796
"10" 29 "Sample2" 1018
"11" 30 "Sample1" 1009
"12" 30 "Sample2" 1366
"13" 31 "Sample1" 1225
"14" 31 "Sample2" 1653
"15" 32 "Sample1" 1579
"16" 32 "Sample2" 2135
"17" 33 "Sample1" 1817
"18" 33 "Sample2" 2678
"19" 34 "Sample1" 2322
"20" 34 "Sample2" 3424
"21" 35 "Sample1" 2861
"22" 35 "Sample2" 4167
"23" 36 "Sample1" 3460
"24" 36 "Sample2" 5220
"25" 37 "Sample1" 4161
"26" 37 "Sample2" 6088
"27" 38 "Sample1" 4930
"28" 38 "Sample2" 7609
"29" 39 "Sample1" 5910
"30" 39 "Sample2" 9308
"31" 40 "Sample1" 7007
"32" 40 "Sample2" 11062
"33" 41 "Sample1" 8339
"34" 41 "Sample2" 13341
"35" 42 "Sample1" 10105
"36" 42 "Sample2" 16288
"37" 43 "Sample1" 11720
"38" 43 "Sample2" 19023
"39" 44 "Sample1" 13312
"40" 44 "Sample2" 22152
"41" 45 "Sample1" 15430
"42" 45 "Sample2" 25360
"43" 46 "Sample1" 17515
"44" 46 "Sample2" 28789
"45" 47 "Sample1" 20022
"46" 47 "Sample2" 32985
"47" 48 "Sample1" 22825
"48" 48 "Sample2" 38019
"49" 49 "Sample1" 26150
"50" 49 "Sample2" 42972
"51" 50 "Sample1" 28998
"52" 50 "Sample2" 49293
"53" 51 "Sample1" 33236
"54" 51 "Sample2" 56239
"55" 52 "Sample1" 37302
"56" 52 "Sample2" 62791
"57" 53 "Sample1" 41212
"58" 53 "Sample2" 69865
"59" 54 "Sample1" 45683
"60" 54 "Sample2" 76922
"61" 55 "Sample1" 49844
"62" 55 "Sample2" 82906
"63" 56 "Sample1" 53833
"64" 56 "Sample2" 90092
"65" 57 "Sample1" 58486
"66" 57 "Sample2" 98194
"67" 58 "Sample1" 64650
"68" 58 "Sample2" 105602
"69" 59 "Sample1" 69514
"70" 59 "Sample2" 114319
"71" 60 "Sample1" 76202
"72" 60 "Sample2" 123807
"73" 61 "Sample1" 83273
"74" 61 "Sample2" 133289
"75" 62 "Sample1" 90783
"76" 62 "Sample2" 146387
"77" 63 "Sample1" 98805
"78" 63 "Sample2" 159302
"79" 64 "Sample1" 107247
"80" 64 "Sample2" 168856
"81" 65 "Sample1" 115200
"82" 65 "Sample2" 180408
"83" 66 "Sample1" 123359
"84" 66 "Sample2" 191083
"85" 67 "Sample1" 134258
"86" 67 "Sample2" 200708
"87" 68 "Sample1" 142179
"88" 68 "Sample2" 212830
"89" 69 "Sample1" 155549
"90" 69 "Sample2" 226788
"91" 70 "Sample1" 165425
"92" 70 "Sample2" 239063
"93" 71 "Sample1" 178681
"94" 71 "Sample2" 255338
"95" 72 "Sample1" 194293
"96" 72 "Sample2" 272725
"97" 73 "Sample1" 211417
"98" 73 "Sample2" 288512
"99" 74 "Sample1" 225639
"100" 74 "Sample2" 307752
"101" 75 "Sample1" 242618
"102" 75 "Sample2" 324919
"103" 76 "Sample1" 259391
"104" 76 "Sample2" 338429
"105" 77 "Sample1" 274041
"106" 77 "Sample2" 354866
"107" 78 "Sample1" 295365
"108" 78 "Sample2" 373456
"109" 79 "Sample1" 315802
"110" 79 "Sample2" 389945
"111" 80 "Sample1" 337826
"112" 80 "Sample2" 411033
"113" 81 "Sample1" 363906
"114" 81 "Sample2" 437069
"115" 82 "Sample1" 388803
"116" 82 "Sample2" 460419
"117" 83 "Sample1" 416728
"118" 83 "Sample2" 487959
"119" 84 "Sample1" 446797
"120" 84 "Sample2" 518295
"121" 85 "Sample1" 474343
"122" 85 "Sample2" 538191
"123" 86 "Sample1" 496615
"124" 86 "Sample2" 558646
"125" 87 "Sample1" 524080
"126" 87 "Sample2" 583043
"127" 88 "Sample1" 549720
"128" 88 "Sample2" 601107
"129" 89 "Sample1" 573864
"130" 89 "Sample2" 623715
"131" 90 "Sample1" 604559
"132" 90 "Sample2" 648209
"133" 91 "Sample1" 632324
"134" 91 "Sample2" 668883
"135" 92 "Sample1" 657924
"136" 92 "Sample2" 698291
"137" 93 "Sample1" 686875
"138" 93 "Sample2" 727313
"139" 94 "Sample1" 712778
"140" 94 "Sample2" 745380
"141" 95 "Sample1" 732786
"142" 95 "Sample2" 766930
"143" 96 "Sample1" 755081
"144" 96 "Sample2" 784337
"145" 97 "Sample1" 772734
"146" 97 "Sample2" 793358
"147" 98 "Sample1" 787861
"148" 98 "Sample2" 804737
"149" 99 "Sample1" 807004
"150" 99 "Sample2" 820443
"151" 100 "Sample1" 823211
"152" 100 "Sample2" 826114
"153" 101 "Sample1" 836362
"154" 101 "Sample2" 839164
"155" 102 "Sample1" 853124
"156" 102 "Sample2" 853144
"157" 103 "Sample1" 865378
"158" 103 "Sample2" 861670
"159" 104 "Sample1" 874578
"160" 104 "Sample2" 875069
"161" 105 "Sample1" 887322
"162" 105 "Sample2" 883285
"163" 106 "Sample1" 888977
"164" 106 "Sample2" 875065
"165" 107 "Sample1" 889700
"166" 107 "Sample2" 877648
"167" 108 "Sample1" 893138
"168" 108 "Sample2" 875875
"169" 109 "Sample1" 887592
"170" 109 "Sample2" 865221
"171" 110 "Sample1" 888684
"172" 110 "Sample2" 863268
"173" 111 "Sample1" 887037
"174" 111 "Sample2" 863513
"175" 112 "Sample1" 887121
"176" 112 "Sample2" 855519
"177" 113 "Sample1" 883085
"178" 113 "Sample2" 854312
"179" 114 "Sample1" 886536
"180" 114 "Sample2" 855445
"181" 115 "Sample1" 875703
"182" 115 "Sample2" 844795
"183" 116 "Sample1" 870749
"184" 116 "Sample2" 837992
"185" 117 "Sample1" 863969
"186" 117 "Sample2" 830157
"187" 118 "Sample1" 849188
"188" 118 "Sample2" 809476
"189" 119 "Sample1" 840690
"190" 119 "Sample2" 798965
"191" 120 "Sample1" 834472
"192" 120 "Sample2" 791200
"193" 121 "Sample1" 822558
"194" 121 "Sample2" 776191
"195" 122 "Sample1" 811864
"196" 122 "Sample2" 765913
"197" 123 "Sample1" 805300
"198" 123 "Sample2" 761200
"199" 124 "Sample1" 794602
"200" 124 "Sample2" 747483
"201" 125 "Sample1" 784650
"202" 125 "Sample2" 739639
"203" 126 "Sample1" 775107
"204" 126 "Sample2" 730152
"205" 127 "Sample1" 763282
"206" 127 "Sample2" 713285
"207" 128 "Sample1" 748652
"208" 128 "Sample2" 699446
"209" 129 "Sample1" 736764
"210" 129 "Sample2" 686165
"211" 130 "Sample1" 723134
"212" 130 "Sample2" 671281
"213" 131 "Sample1" 710686
"214" 131 "Sample2" 658314
"215" 132 "Sample1" 701834
"216" 132 "Sample2" 650564
"217" 133 "Sample1" 689531
"218" 133 "Sample2" 635331
"219" 134 "Sample1" 678070
"220" 134 "Sample2" 628475
"221" 135 "Sample1" 668757
"222" 135 "Sample2" 617210
"223" 136 "Sample1" 654506
"224" 136 "Sample2" 601336
"225" 137 "Sample1" 644720
"226" 137 "Sample2" 592700
"227" 138 "Sample1" 630305
"228" 138 "Sample2" 579707
"229" 139 "Sample1" 615571
"230" 139 "Sample2" 562022
"231" 140 "Sample1" 603896
"232" 140 "Sample2" 550812
"233" 141 "Sample1" 597816
"234" 141 "Sample2" 540016
"235" 142 "Sample1" 584121
"236" 142 "Sample2" 526066
"237" 143 "Sample1" 571689
"238" 143 "Sample2" 516774
"239" 144 "Sample1" 564142
"240" 144 "Sample2" 509033
"241" 145 "Sample1" 552281
"242" 145 "Sample2" 496389
"243" 146 "Sample1" 543184
"244" 146 "Sample2" 487283
"245" 147 "Sample1" 532683
"246" 147 "Sample2" 477035
"247" 148 "Sample1" 520019
"248" 148 "Sample2" 461698
"249" 149 "Sample1" 510484
"250" 149 "Sample2" 451095
"251" 150 "Sample1" 507832
"252" 150 "Sample2" 446390
"253" 151 "Sample1" 495343
"254" 151 "Sample2" 431015
"255" 152 "Sample1" 484169
"256" 152 "Sample2" 417551
"257" 153 "Sample1" 471422
"258" 153 "Sample2" 404674
"259" 154 "Sample1" 460944
"260" 154 "Sample2" 390280
"261" 155 "Sample1" 448227
"262" 155 "Sample2" 380134
"263" 156 "Sample1" 437738
"264" 156 "Sample2" 369703
"265" 157 "Sample1" 428359
"266" 157 "Sample2" 354902
"267" 158 "Sample1" 417445
"268" 158 "Sample2" 343343
"269" 159 "Sample1" 408065
"270" 159 "Sample2" 331193
"271" 160 "Sample1" 395244
"272" 160 "Sample2" 315775
"273" 161 "Sample1" 386510
"274" 161 "Sample2" 305048
"275" 162 "Sample1" 377028
"276" 162 "Sample2" 292336
"277" 163 "Sample1" 366563
"278" 163 "Sample2" 277876
"279" 164 "Sample1" 357659
"280" 164 "Sample2" 267642
"281" 165 "Sample1" 348242
"282" 165 "Sample2" 255871
"283" 166 "Sample1" 337553
"284" 166 "Sample2" 242636
"285" 167 "Sample1" 327824
"286" 167 "Sample2" 231054
"287" 168 "Sample1" 317584
"288" 168 "Sample2" 219066
"289" 169 "Sample1" 306454
"290" 169 "Sample2" 204829
"291" 170 "Sample1" 295255
"292" 170 "Sample2" 192423
"293" 171 "Sample1" 284770
"294" 171 "Sample2" 180440
"295" 172 "Sample1" 274444
"296" 172 "Sample2" 166529
"297" 173 "Sample1" 264469
"298" 173 "Sample2" 155605
"299" 174 "Sample1" 253923
"300" 174 "Sample2" 145166
"301" 175 "Sample1" 243462
"302" 175 "Sample2" 134932
"303" 176 "Sample1" 233021
"304" 176 "Sample2" 124474
"305" 177 "Sample1" 223360
"306" 177 "Sample2" 114975
"307" 178 "Sample1" 213341
"308" 178 "Sample2" 105267
"309" 179 "Sample1" 202191
"310" 179 "Sample2" 95632
"311" 180 "Sample1" 193325
"312" 180 "Sample2" 87027
"313" 181 "Sample1" 182022
"314" 181 "Sample2" 78802
"315" 182 "Sample1" 173011
"316" 182 "Sample2" 71116
"317" 183 "Sample1" 164207
"318" 183 "Sample2" 64758
"319" 184 "Sample1" 154985
"320" 184 "Sample2" 57082
"321" 185 "Sample1" 147070
"322" 185 "Sample2" 51728
"323" 186 "Sample1" 138927
"324" 186 "Sample2" 46295
"325" 187 "Sample1" 130203
"326" 187 "Sample2" 40908
"327" 188 "Sample1" 122935
"328" 188 "Sample2" 36275
"329" 189 "Sample1" 115565
"330" 189 "Sample2" 31825
"331" 190 "Sample1" 107969
"332" 190 "Sample2" 28028
"333" 191 "Sample1" 100995
"334" 191 "Sample2" 24447
"335" 192 "Sample1" 94510
"336" 192 "Sample2" 21125
"337" 193 "Sample1" 88365
"338" 193 "Sample2" 18489
"339" 194 "Sample1" 82247
"340" 194 "Sample2" 15878
"341" 195 "Sample1" 77315
"342" 195 "Sample2" 13881
"343" 196 "Sample1" 71657
"344" 196 "Sample2" 12183
"345" 197 "Sample1" 66920
"346" 197 "Sample2" 10305
"347" 198 "Sample1" 61298
"348" 198 "Sample2" 9096
"349" 199 "Sample1" 57765
"350" 199 "Sample2" 7759
"351" 200 "Sample1" 53811
"352" 200 "Sample2" 6645
"353" 201 "Sample1" 48906
"354" 201 "Sample2" 5532
"355" 202 "Sample1" 44847
"356" 202 "Sample2" 4642
"357" 203 "Sample1" 41546
"358" 203 "Sample2" 3997
"359" 204 "Sample1" 38125
"360" 204 "Sample2" 3360
"361" 205 "Sample1" 35039
"362" 205 "Sample2" 2858
"363" 206 "Sample1" 31902
"364" 206 "Sample2" 2447
"365" 207 "Sample1" 29076
"366" 207 "Sample2" 2185
"367" 208 "Sample1" 26271
"368" 208 "Sample2" 1811
"369" 209 "Sample1" 24037
"370" 209 "Sample2" 1494
"371" 210 "Sample1" 21487
"372" 210 "Sample2" 1313
"373" 211 "Sample1" 19187
"374" 211 "Sample2" 1028
"375" 212 "Sample1" 17317
"376" 212 "Sample2" 954
"377" 213 "Sample1" 15741
"378" 213 "Sample2" 819
"379" 214 "Sample1" 14000
"380" 214 "Sample2" 693
"381" 215 "Sample1" 12494
"382" 215 "Sample2" 614
"383" 216 "Sample1" 10963
"384" 216 "Sample2" 556
"385" 217 "Sample1" 9928
"386" 217 "Sample2" 549
"387" 218 "Sample1" 8652
"388" 218 "Sample2" 421
"389" 219 "Sample1" 7604
"390" 219 "Sample2" 404
"391" 220 "Sample1" 6515
"392" 220 "Sample2" 394
"393" 221 "Sample1" 5735
"394" 221 "Sample2" 325
"395" 222 "Sample1" 4783
"396" 222 "Sample2" 301
"397" 223 "Sample1" 4108
"398" 223 "Sample2" 308
"399" 224 "Sample1" 3611
"400" 224 "Sample2" 290
"401" 225 "Sample1" 3224
"402" 225 "Sample2" 307
"403" 226 "Sample1" 2676
"404" 226 "Sample2" 304
"405" 227 "Sample1" 2261
"406" 227 "Sample2" 273
"407" 228 "Sample1" 1936
"408" 228 "Sample2" 303
"409" 229 "Sample1" 1664
"410" 229 "Sample2" 302
"411" 230 "Sample1" 1529
"412" 230 "Sample2" 291
"413" 231 "Sample1" 1180
"414" 231 "Sample2" 313
"415" 232 "Sample1" 975
"416" 232 "Sample2" 288
"417" 233 "Sample1" 897
"418" 233 "Sample2" 274
"419" 234 "Sample1" 812
"420" 234 "Sample2" 281
"421" 235 "Sample1" 712
"422" 235 "Sample2" 318
"423" 236 "Sample1" 587
"424" 236 "Sample2" 288
"425" 237 "Sample1" 555
"426" 237 "Sample2" 299
"427" 238 "Sample1" 524
"428" 238 "Sample2" 314
"429" 239 "Sample1" 443
"430" 239 "Sample2" 336
"431" 240 "Sample1" 411
"432" 240 "Sample2" 354
"433" 241 "Sample1" 382
"434" 241 "Sample2" 339
"435" 242 "Sample1" 344
"436" 242 "Sample2" 342
"437" 243 "Sample1" 340
"438" 243 "Sample2" 366
"439" 244 "Sample1" 362
"440" 244 "Sample2" 353
"441" 245 "Sample1" 325
"442" 245 "Sample2" 373
"443" 246 "Sample1" 290
"444" 246 "Sample2" 340
"445" 247 "Sample1" 319
"446" 247 "Sample2" 345
"447" 248 "Sample1" 293
"448" 248 "Sample2" 341
"449" 249 "Sample1" 304
"450" 249 "Sample2" 365
"451" 250 "Sample1" 354
"452" 250 "Sample2" 381
"453" 251 "Sample1" 296
"454" 251 "Sample2" 366
"455" 252 "Sample1" 315
"456" 252 "Sample2" 415
"457" 253 "Sample1" 315
"458" 253 "Sample2" 399
"459" 254 "Sample1" 331
"460" 254 "Sample2" 377
"461" 255 "Sample1" 321
"462" 255 "Sample2" 410
"463" 256 "Sample1" 326
"464" 256 "Sample2" 421
"465" 257 "Sample1" 349
"466" 257 "Sample2" 433
"467" 258 "Sample1" 370
"468" 258 "Sample2" 429
"469" 259 "Sample1" 318
"470" 259 "Sample2" 439
"471" 260 "Sample1" 344
"472" 260 "Sample2" 433
"473" 261 "Sample1" 346
"474" 261 "Sample2" 400
"475" 262 "Sample1" 341
"476" 262 "Sample2" 419
"477" 263 "Sample1" 345
"478" 263 "Sample2" 499
"479" 264 "Sample1" 348
"480" 264 "Sample2" 418
"481" 265 "Sample1" 371
"482" 265 "Sample2" 454
"483" 266 "Sample1" 405
"484" 266 "Sample2" 441
"485" 267 "Sample1" 380
"486" 267 "Sample2" 476
"487" 268 "Sample1" 365
"488" 268 "Sample2" 464
"489" 269 "Sample1" 366
"490" 269 "Sample2" 461
"491" 270 "Sample1" 351
"492" 270 "Sample2" 483
"493" 271 "Sample1" 431
"494" 271 "Sample2" 481
"495" 272 "Sample1" 340
"496" 272 "Sample2" 503
"497" 273 "Sample1" 405
"498" 273 "Sample2" 493
"499" 274 "Sample1" 370
"500" 274 "Sample2" 529
"501" 275 "Sample1" 374
"502" 275 "Sample2" 469
"503" 276 "Sample1" 385
"504" 276 "Sample2" 503
"505" 277 "Sample1" 361
"506" 277 "Sample2" 496
"507" 278 "Sample1" 410
"508" 278 "Sample2" 491
"509" 279 "Sample1" 410
"510" 279 "Sample2" 521
"511" 280 "Sample1" 368
"512" 280 "Sample2" 514
"513" 281 "Sample1" 397
"514" 281 "Sample2" 501
"515" 282 "Sample1" 379
"516" 282 "Sample2" 546
"517" 283 "Sample1" 384
"518" 283 "Sample2" 502
"519" 284 "Sample1" 389
"520" 284 "Sample2" 513
"521" 285 "Sample1" 419
"522" 285 "Sample2" 571
"523" 286 "Sample1" 403
"524" 286 "Sample2" 519
"525" 287 "Sample1" 423
"526" 287 "Sample2" 536
"527" 288 "Sample1" 480
"528" 288 "Sample2" 653
"529" 289 "Sample1" 478
"530" 289 "Sample2" 564

Here's my attempt.
library(tidyverse)
library(S4Vectors)

# This would probably be read in using `read_tsv()` or something
rawData <- structure(
    list(
        Readlength = c(25, 26, 27, 28, 29, 30, 31, 32,
                       33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48,
                       49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64,
                       65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
                       81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96,
                       97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109,
                       110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122,
                       123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135,
                       136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148,
                       149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161,
                       162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174,
                       175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187,
                       188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200,
                       201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213,
                       214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226,
                       227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239,
                       240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252,
                       253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265,
                       266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278,
                       279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289),
        Sample1 = c(457, 603, 599, 686, 796, 1009, 1225, 1579, 1817, 2322, 2861, 3460,
                    4161, 4930, 5910, 7007, 8339, 10105, 11720, 13312, 15430, 17515,
                    20022, 22825, 26150, 28998, 33236, 37302, 41212, 45683, 49844,
                    53833, 58486, 64650, 69514, 76202, 83273, 90783, 98805, 107247,
                    115200, 123359, 134258, 142179, 155549, 165425, 178681, 194293,
                    211417, 225639, 242618, 259391, 274041, 295365, 315802, 337826,
                    363906, 388803, 416728, 446797, 474343, 496615, 524080, 549720,
                    573864, 604559, 632324, 657924, 686875, 712778, 732786, 755081,
                    772734, 787861, 807004, 823211, 836362, 853124, 865378, 874578,
                    887322, 888977, 889700, 893138, 887592, 888684, 887037, 887121,
                    883085, 886536, 875703, 870749, 863969, 849188, 840690, 834472,
                    822558, 811864, 805300, 794602, 784650, 775107, 763282, 748652,
                    736764, 723134, 710686, 701834, 689531, 678070, 668757, 654506,
                    644720, 630305, 615571, 603896, 597816, 584121, 571689, 564142,
                    552281, 543184, 532683, 520019, 510484, 507832, 495343, 484169,
                    471422, 460944, 448227, 437738, 428359, 417445, 408065, 395244,
                    386510, 377028, 366563, 357659, 348242, 337553, 327824, 317584,
                    306454, 295255, 284770, 274444, 264469, 253923, 243462, 233021,
                    223360, 213341, 202191, 193325, 182022, 173011, 164207, 154985,
                    147070, 138927, 130203, 122935, 115565, 107969, 100995, 94510,
                    88365, 82247, 77315, 71657, 66920, 61298, 57765, 53811, 48906,
                    44847, 41546, 38125, 35039, 31902, 29076, 26271, 24037, 21487,
                    19187, 17317, 15741, 14000, 12494, 10963, 9928, 8652, 7604, 6515,
                    5735, 4783, 4108, 3611, 3224, 2676, 2261, 1936, 1664, 1529, 1180,
                    975, 897, 812, 712, 587, 555, 524, 443, 411, 382, 344, 340, 362,
                    325, 290, 319, 293, 304, 354, 296, 315, 315, 331, 321, 326, 349,
                    370, 318, 344, 346, 341, 345, 348, 371, 405, 380, 365, 366, 351,
                    431, 340, 405, 370, 374, 385, 361, 410, 410, 368, 397, 379, 384,
                    389, 419, 403, 423, 480, 478),
        Sample2 = c(580, 648, 715, 914,
                    1018, 1366, 1653, 2135, 2678, 3424, 4167, 5220, 6088, 7609, 9308,
                    11062, 13341, 16288, 19023, 22152, 25360, 28789, 32985, 38019,
                    42972, 49293, 56239, 62791, 69865, 76922, 82906, 90092, 98194,
                    105602, 114319, 123807, 133289, 146387, 159302, 168856, 180408,
                    191083, 200708, 212830, 226788, 239063, 255338, 272725, 288512,
                    307752, 324919, 338429, 354866, 373456, 389945, 411033, 437069,
                    460419, 487959, 518295, 538191, 558646, 583043, 601107, 623715,
                    648209, 668883, 698291, 727313, 745380, 766930, 784337, 793358,
                    804737, 820443, 826114, 839164, 853144, 861670, 875069, 883285,
                    875065, 877648, 875875, 865221, 863268, 863513, 855519, 854312,
                    855445, 844795, 837992, 830157, 809476, 798965, 791200, 776191,
                    765913, 761200, 747483, 739639, 730152, 713285, 699446, 686165,
                    671281, 658314, 650564, 635331, 628475, 617210, 601336, 592700,
                    579707, 562022, 550812, 540016, 526066, 516774, 509033, 496389,
                    487283, 477035, 461698, 451095, 446390, 431015, 417551, 404674,
                    390280, 380134, 369703, 354902, 343343, 331193, 315775, 305048,
                    292336, 277876, 267642, 255871, 242636, 231054, 219066, 204829,
                    192423, 180440, 166529, 155605, 145166, 134932, 124474, 114975,
                    105267, 95632, 87027, 78802, 71116, 64758, 57082, 51728, 46295,
                    40908, 36275, 31825, 28028, 24447, 21125, 18489, 15878, 13881,
                    12183, 10305, 9096, 7759, 6645, 5532, 4642, 3997, 3360, 2858,
                    2447, 2185, 1811, 1494, 1313, 1028, 954, 819, 693, 614, 556,
                    549, 421, 404, 394, 325, 301, 308, 290, 307, 304, 273, 303, 302,
                    291, 313, 288, 274, 281, 318, 288, 299, 314, 336, 354, 339, 342,
                    366, 353, 373, 340, 345, 341, 365, 381, 366, 415, 399, 377, 410,
                    421, 433, 429, 439, 433, 400, 419, 499, 418, 454, 441, 476, 464,
                    461, 483, 481, 503, 493, 529, 469, 503, 496, 491, 521, 514, 501,
                    546, 502, 513, 571, 519, 536, 653, 564)),
    class = c("spec_tbl_df", "tbl_df", "tbl", "data.frame"),
    row.names = c(NA, -265L),
    spec = structure(
        list(
            cols = list(
                Readlength = structure(
                    list(), class = c("collector_double", "collector")
                    ),
                Sample1 = structure(
                    list(), class = c("collector_double", "collector")
                ),
                Sample2 = structure(
                    list(), class = c("collector_double", "collector"))
                ),
            default = structure(
                list(), class = c("collector_guess", "collector")
            ),
            skip = 1),
        class = "col_spec")
    )

rlData <- rawData %>%
    pivot_longer( # Put all the counts in a single column
        cols = contains("Sample"),
        names_to = "Sample",
        values_to = "Count"
        ) %>%
    split(f = .$Sample) %>% # Separate into each sample
lapply(function(x){
    # Now create a run-length encoded vector
    # This encodes millions of repetitive values into 2 numbers
    # 1) the value & 2) how many
    Rle(values = x$Readlength, lengths = x$Count)
        # They can be a bit difficult to plot though
    })

# I would create the quantiles for the boxplot separately now
q <- c(min = 0, q25 = 0.25, median = 0.5, q75 = 0.75, max = 1)
# Collect the stats
bxp <- rlData %>%
    lapply(quantile, probs = q) %>% # for every Rle, get the stats
    as_tibble() %>%
    mutate(stat = names(q)) %>% # This just collates the output
    pivot_longer(
        # Make it long
        cols = contains("Sample"),
        names_to = "Sample",
        values_to = "value"
    ) %>%
    pivot_wider(
        # Now spread it out for each stat/column
        names_from = stat,
        values_from = value
    ) %>%
    mutate(
        # set a few values for easier plotting
        IQR = q75 - q25,
        Sample = as.factor(Sample),
        x = as.integer(Sample)
    )

w <- 0.4 # Box width
bxp %>%
    mutate(
        # Set the width of the boxes
        xmin = x - w,
        xmax = x + w,
        # Set the values for the whiskers.
        # These usually extend 1.5*IQR.
        # Any points beyond are considered outliers
        boxMax = case_when(
            max > q75 + 1.5*IQR ~ q75 + 1.5*IQR,
            max > q75 + 1.5*IQR ~ max
        ),
        boxMin = case_when(
            min < q25 - 1.5*IQR ~ q25 - 1.5*IQR,
            min > q25 - 1.5*IQR ~ min
        )
    ) %>%
    ggplot(aes(x = x, y = median)) +
    geom_rect(
        # Draw the box
        aes(xmin = xmin, xmax = xmax, ymin = q25, ymax = q75),
        fill = "grey80",
        colour = "black"
    ) +
    geom_segment(
        # Add the median
        aes(x = xmin, xend = xmax, yend = median),
        colour = "black"
    ) +
    geom_segment(
        # The top whiskers
        aes(x = x, xend = x, y = q75, yend = boxMax)
    ) +
    geom_segment(
        # The bottom whiskers
        aes(x = x, xend = x, y = q25, yend = boxMin)
    ) +
    geom_point(aes(y = max)) + # Add some symbolic outliers if you want
    geom_point(aes(y = min)) +
    scale_x_continuous(
        # Tidy up the x-axis ticks
        breaks = bxp$x,
    labels = bxp$Sample
    ) +
    labs(
        # Tidy up the labels
        x = "Library",
        y = "Read Length"
    ) +
    theme_bw() # Get rid of the stupid grey background

That gives me this figure

Hope that solves it for you. I probably went the long way around though :)

Just what I wanted, thank you!

My first recommendation would be to use (i.e. ask for) the indidividual-level data, which is what ggplot2 expects.
Binned Data
If that's not possible, binned data can be shoehorned into a ggplot2 boxplot by forcing it back into individual values. My example here adds new unbinned rows to the table with a count of -1, then removes rows with a positive count:
library(tidyverse)
library(ggplot2)

read_table2("11196_input.txt") %>%
    pivot_longer(cols=c("Sample1", "Sample2"),
                 names_to="Sample", values_to="Count") %>%
    bind_rows(tibble(Length=rep(.$Length, .$Count),
                     Sample=rep(.$Sample, .$Count),
                     Count=rep(.$Length * 0 - 1, .$Count))) %>%
    filter(Count < 0) %>%
    select(-Count) -> data.tbl

data.tbl %>%
    ggplot() +
    aes(x = Sample, y = Length) +
    geom_boxplot()

It would be possible to make the boxplot look more normal by modifying the rep functions to create values based on a distribution (e.g. using runif or rnorm), but that gets quite tricky because (for example) binned normal data would have different probability distributions for each bin.
Individual Values
If, on the other hand, individual values are available in long format, creating a boxplot becomes much easier because the data fits the format expected by ggplot2:
library(tidyverse)
library(ggplot2)

read.table("11196_input2.txt") -> data.tbl

data.tbl %>%
    ggplot() +
    aes(x = name, y = Readlength) +
    geom_boxplot()

I don't think that the individual example worked -- the box plots should be different as both of the samples have different underlying data.

Ah, yeah. That's because the data isn't individual-level (i.e. one line per observation), and I didn't realise that the heading change from "count" to "value" still represented the same thing. I'm seeing 130M data points... that's quite a lot for ggplot2 to process. It'd probably be better to use stat_boxplot rather than geom_boxplot, and specify the boxplot points rather than the values.

The plot you pasted in your question has two layers of information: i) fragment (or read) length, ii) sample.
To be able to draw the same plot with your data using ggplot2, you will need your data to be in the "long format" (your data as you presented now is in the "wide format"):
library(readr)
library(tidyr)
library(ggplot2)

my_data <- "
Length Sample1 Sample2
20 10 5
30 15 223
40 16 4923
50 293 239
60 3290 23
70 248 6
80 23 1
90 21 0
100 5 0
"

my_data <- read_delim(my_data, delim = " ")

my_data_long <- pivot_longer(my_data, cols = c("Sample1", "Sample2"))

> my_data_long
# A tibble: 18 x 3
   Length name    value
    <dbl> <chr>   <dbl>
 1     20 Sample1    10
 2     20 Sample2     5
 3     30 Sample1    15
 4     30 Sample2   223
 5     40 Sample1    16
 6     40 Sample2  4923
 7     50 Sample1   293
 8     50 Sample2   239
 9     60 Sample1  3290
10     60 Sample2    23
11     70 Sample1   248
12     70 Sample2     6
13     80 Sample1    23
14     80 Sample2     1
15     90 Sample1    21
16     90 Sample2     0
17    100 Sample1     5
18    100 Sample2     0

Moreover, boxplots are meaningful when you would like to visualize distributions, in your case you don't have all of your data points but only a summary of those, counts of read lengths. The best I could think of was using a third aesthetic to reflect the counts (after removing read lengths that were not observed):
# remove read lengths that are not observed
my_data_long <- my_data_long[my_data_long$value != 0,]

ggplot(my_data_long,
       aes(x = name, y = Length)) +
  geom_boxplot() +
  geom_point(aes(size = value))

EDIT: In a boxplot, the quantiles reflect the number of data points but the way I draw the boxplot (using Length and name columns) do not make use of such information. A better way of visualising the data would be the following (and is not that different from what @gringer suggested):
ggplot(my_data_long,
       aes(x = name, y = Length)) +
  geom_point(aes(size = value))


text:
CLI access to genomes/proteomes by BioSample ID
This question was also asked on Biostars
I'd like to download multiple genome assemblies or proteomes using a set of BioSample IDs from NCBI.
I'm able to find the assemblies belonging to the BioSample IDS in a browser (in the search field of https://www.ncbi.nlm.nih.gov/), but couldn't find a commandline interface solution (CLI).
E.g. for BioSample SAMN09405588 the assembly id is PDT000806148.1, and from https://www.ncbi.nlm.nih.gov/assembly/GCA_014136285.1/ I can download the proteome: GCA_014136285.1_PDT000806148.1_protein.faa.gz
Thank you for your help!

I think this question is not a duplicate of https://bioinformatics.stackexchange.com/questions/15794/using-entrez-efetch-to-retrive-fasta-file-from-any-ncbi-database - there a gene id is given, here a biosample id.

Answer from @m, converted from comment:
Biostars are recommending esearch ... its not the approach I'd have taken but presume will do the job. Its available via SRA, so I'd have gone for the SRA toolkit. Is SRA more appropriate in context? (I liked Biostars BeautifulSoup solution)
[Please update or improve this answer if you can]

source: Stack Exchange

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

sort: votes


text:
Finding insertion in a plasmid
This question was also asked on Biostars
I received a set of Nanopore long reads of a ~3kbp long known plasmid with an unknown insert between 1kbp and 7kbp. Could you suggest a pipeline to assemble the reads and find the insertion?

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.

When I've needed to do this, I use canu to assemble the nanopore reads (choosing the highest-quality reads that are at least half the length of the target), then map the assembled consensus back to the reference sequence.
The Canu assembly part of that is described here:
https://www.protocols.io/view/plasmid-sequence-analysis-from-long-reads-36wgq4n5yvk5?step=12
There are numerous different ways to map the consensus assembled sequence to the reference sequence. BLASTn in "align 2 or more sequences" mode could be used, for example. Another option, because the sequences are relatively small and you're just aligning 2, is needle from the EMBOSS suite.

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


text:
loading in .txt files with different spacing in column names into R
Trying to read in this file but saved as a .txt file: http://www.informatics.jax.org/downloads/reports/HOM_MouseHumanSequence.rpt
gene_conversion = read.table("HOM_MouseHumanSequence.rpt.txt",header=T, fill=T)
but the column names are splitting weirdly.
> head(gene_conversion)
        DB  Class        Key Common Organism       Name        NCBI Taxon    ID                      Symbol
1 38710983 mouse, laboratory  10090     Gdnf      14573  MGI:107430 Chr15   3.8                          cM
2 38710983  human       9606   GDNF     2668  HGNC:4232 OMIM:600837  Chr5 p13.2   Chr5:37812677-37840044(-)
3 38710984 mouse, laboratory  10090     Rrm1      20133   MGI:98180  Chr7 54.72                          cM
4 38710984  human       9606   RRM1     6240 HGNC:10451 OMIM:180410 Chr11 p15.4    Chr11:4094685-4138925(+)
5 38710985 mouse, laboratory  10090   Hoxd12      15432   MGI:96204  Chr2 44.13                          cM
6 38710985  human       9606 HOXD12     3238  HGNC:5135 OMIM:142988  Chr2 q31.1 Chr2:176099802-176100760(+)
                                                            EntrezGene
1                                             Chr15:7840327-7867056(+)
2 NM_199231,NM_001278098,NM_199234,NM_001190469,NM_001190468,NM_000514
3                                          Chr7:102090902-102118978(+)
4                     NM_001033,NM_001318065,NM_001330193,NM_001318064
5                                            Chr2:74505357-74508049(+)
6                                                            NM_021193
                                                                                ID.1
1                                     NM_001301357,NM_010275,NM_001301333,NM_001301332
2 NP_001265027,NP_000505,NP_001177397,NP_001177398,XP_016864826,XP_011512332,NP_954701
3                                                                            NM_009103
4                                     NP_001024,NP_001317122,NP_001304994,NP_001304993
5                                                                            NM_008274
6                                                                            NP_067016
                                             Mouse    MGI ID.2 HGNC ID.3 OMIM Gene ID.4 Genetic Location Genomic
1 NP_001288286,NP_001288262,NP_034405,NP_001288261 P48540   NA   NA   NA   NA   NA   NA      NA       NA      NA
2                                           P39905          NA   NA   NA   NA   NA   NA      NA       NA      NA
3                                        NP_033129 P07742   NA   NA   NA   NA   NA   NA      NA       NA      NA
4                                           P23921          NA   NA   NA   NA   NA   NA      NA       NA      NA
5                                        NP_032300 P23812   NA   NA   NA   NA   NA   NA      NA       NA      NA
6                                           P35452          NA   NA   NA   NA   NA   NA      NA       NA      NA
  Coordinates X.mouse. X. human. X..1 Nucleotide RefSeq IDs Protein RefSeq.1 IDs.1 SWISS_PROT IDs.2
1          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA
2          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA
3          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA
4          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA
5          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA
6          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA

I see spaces in the headings, which suggests this might be a tab-separated file, rather than a space-separated file. If that's the case, you should be using sep="\t" in your read.table command, or alternatively using read.delim:
gene_conversion <-
  read.table("HOM_MouseHumanSequence.rpt.txt", sep="\t", header=T, fill=T)

or:
gene_conversion <- read.delim("HOM_MouseHumanSequence.rpt.txt")

On a related note, how can I subset based on just a column?

Can you please ask a separate question for that?


text:
How to calculate coverage depth of genes in an annotated sequence?
I have a metagenomic dataset across multiple samples that has been grouped and binned (we did this due to poor sequence quallity).
It is several samples taken from termite guts of differing castes, the samples have been grouped to make binning easier due to poor sequence suitability. They have been binned via ASVs, assembling reads from all samples to generate the putative bacterial species. So I have fasta files of putative species that have been annotated via prodigal. I want to now calculate the abundance of each bin in a given sample. I think the best way to do this would be at the gene level by looking at the number of reads that map to a gene in any given sample. I am just unsure how to do this.
I have annotated the bins and now want to calculate the abundance of each bin in a sample but looking at the number of reads from a sample that map to each gene in a bin.
Is this a viable method? If so how would I do it and if not what else could I do?

The best way to approach it would be to take trimmed reads map them to the bins (i.e. with bowtie2) and estimate the abundance of each bin with samtools by getting the depth of coverage:
https://www.metagenomics.wiki/tools/samtools/breadth-of-coverage
or you can use QUANT_BIN module from the METAWRAP to do it for you. Follow this for the METAWRAP approach:
https://github.com/bxlab/metaWRAP/issues/40


text:
RNA_Seq data aligned used uniquely or multi mapped reads impact on result interpretation
I have some transcriptomic (Whole) sequencing data that I should analyse. I would like to do raw data alignment to a reference genome taking into account the multi mapped reads and uniquely mapped reads. I would like to know the consequences on my transcriptomic analysis results : 1. if I should only focus on uniqueMapping reads (mapped to one loci) and neglect MultiMapping read (mapped to several loci or variants)! 2. or take them all together
So my questions are :
1- Considering unique and multi-aligned reads, wouldn't this be a source of bias in my results interpretation?
2- By hiding the multi-mapped reads, could I lose information on the expression of important genes?
thank you in advance

There are several methods/tools that take into account multi-mappers and deal with them in reasonable ways (e.g. pseudo-alignment methods like salmon and more traditional alignment-based methods like RSEM). If you completely ignore multi-mappers, then yes, you will be losing information, which may or may not be valuable to you.
Both of those methods utilize the expectation-maximization algorithm (or a variant thereof) to determine the relative contributions of multi-mapping reads to the various mapping locations. With salmon, you can forego alignment completely and end up with very accurate gene quantifications in hand in like 20 minutes if you want.

I followed your advice and I didn't have a big difference between the results of considering single mapped and multimapped reads. Thanks you Jared Andrews


text:
What is the significance of enrichment factor regarding machine learning methods?
I'm fairly new to the field of bioinformatics and ran into a question while reading a paper I found on bioRxiv. The overall setting of the paper is using multi-task deep neural networks for kinase classification. The paper that I'm referring to is Kinome-Wide Activity Classification of Small Molecules by Deep Learning (Allen et al., 2019) for anyone who's curious.
I'll get straight into the specific figure that I'm having trouble understanding.

I know that questions like these are a bit discouraged on SE, but I'm having trouble on how to start regarding the interpretation of this figure.
Based on my understanding, the point seems to be that multitaks DNN's have the highest overall EF/EFmax ratio amounts, but what does this ratio signify?
I did some research and according to Wikipedia an "enrichment factor" in terms of bioinformatics is a concept "to measure the added value of a search tool over another one or over the homogeneous distribution in the genome population." Does this mean that multitask DNN's are "better" than the other models, due to its higher number of EF/EFmax points?
But if we look at the curves on the right, it seems to say that for fewer number of compounds, the ratio for multitask DNN's is lower...

I haven't read the paper, I don't do small molecules at all. The thing that is interesting is the multi-task DNN (deep neural network). Any chance you can explain its architecture please? DNN will possibly outperform the machine learn ing methods architecture and data set size depending. That is the result of this investigation, but more important is how they constructed the DNN, if they have borrowed it from natural language processing - I am not impressed. If they constructed it de novo that is a state-of-the-art approach.

The final question, given this is very likely under review are you reviewing the manuscript? The ratio is their accuracy estimate and is a standard approach.

I did some research and according to Wikipedia an "enrichment factor" in terms of bioinformatics is...

Where possible, it's best to use explanations provided in the original source. This is because definitions could be different, depending on who is writing. Here's some descriptive information from the paper that may be more informative:

The enrichment factor (EF) is another metric we use and is very
  popular for use in evaluation of virtual screening performance. EF
  denotes the quotient of true actives among a subset of predicted
  actives and the overall fraction of actives.... To evaluate the
  maximum achievable enrichment factor (EFmax), the maximum number of
  actives among the percentage of tested compounds is divided by the
  fraction of actives in the entire dataset. To more consistently
  quantify enrichment at 0.1% of all compounds tested, we report the
  ratio of EF/ EFmax. This normalized metric is useful because EF values
  are not directly comparable across different datasets due to the
  maximum possible EF being constrained by the ratio of total to active
  compounds and the evaluated fraction as shown in the equation above.

So for the purpose of this paper, they're saying that this ratio is an appropriately normalised substitute for a plain enrichment factor. High enrichment factor == good; low enrichment factor == bad.
Given the information you've provided, I'd agree that the curves in this figure appear to suggest that the multitask DNN performs worse with smaller numbers of active compounds (where Naive Bayes or Logistic regression would be a better choice). It's hard to tell from the graph (and legend) what the curves represent, and the clipping at EF/EFmax = 0 doesn't help.
The curves have a strange shape, suggesting the number of active compounds is dependent on the EF/EFmax ratio, rather than the other way round (which would make more sense to me). If I consider the curves to be something like a smoothed mean of the active compound number, then any place where there are fewer dots (i.e. low EF/EFmax for Multitask DNN) will have more associated error. For Multitask DNN, I see a large gap between 0.0 and 0.2, with lots of outliers on the 0.0 line (for all methods), so presumably the 0.0 outliers are skewing the curves.
Given that this is a preprint, the generation of the curves isn't specified, and the other figures (especially Figure 3) seem to tell a different story (i.e. that multi-task DNN performs better at all scales), I'd recommend contacting the authors to ask / clarify what the curves represent, or ignoring the curves entirely.

Ok I'll be brief, given the ms is likely under review.
Background in silico drug discovery has been part of deep learning for a long time 2015(?), and tensorflow was deployed. They worked out how to vectorise a small molecule and thereafter the entire field was.wide open to advanced ML.
Ratio The ratio is an accuracy index of the machine learning algorithms. The denominator is the total number of biological observations, the numerator is the predicted number. This 90% would be an accuracy index, i.e. the algorithm is correct 9 out of 10 times.
Result The plot shows, as the data set size increases the accuracy of the prediction increases most notably for the DNN and that conforms to all the theory. DNN/ANN only become powerful for 'big data', standard ML, i.e. random forest, naive Bayes, LR work well on small data sets, are fast and produce a classification which isn't wrapped up in a 3-4 Gig dump, thereafter performance saturates for routine ML methods.
The multi task DNN (MT-DNN) means that more than one loss function is being used, which is unusual in generic deep learning, but has been applied to natural language and drug discovery previously. Why this is being deployed is a long explanation, so I'll leave it there.


text:
Calculating most abundant transcript from RNA-Seq data
vcf2maf uses VEP to annotate variants, and I believe selects the default Ensembl transcript to use for annotation. Sometimes the transcript that VEP selects is not the transcript I'm interested in, usually because the selected transcript is not the most highly expressed transcript in my tissue of interest (skin). vcf2maf allows you to provide a transcript override list so that VEP annotates the variant using the specified transcripts instead.
I have several skin samples sequenced with bulk RNA-Seq. I want to estimate the average abundance for each transcript across all samples and then use these abundances to rank transcripts from most to least abundant. Then I will use the most abundant transcript as the default VEP transcript. I plan to use salmon or kallisto to quantify transcript abundance. Should I use TPM or normalized counts to calculate average expression?
My initial thought is to use normalized counts (generated by DESeq2 from raw counts). Are there any problems with this approach? GTEx displays transcript abundance with average TPM, but I thought TPM was inappropriate to use across samples because it doesn't account for between sample differences.
Update: I forgot to mention I also tried using TPM ranks like @ATpoint describes. I haven't fully compared how this compares to transcripts identified by normalized counts, but the initial genes I checked showed good concordance between methods

Transcript abundance quantification is a tricky topic since a read often could belong to several transcripts, so any "count" is a best guess as to which transcript it actually originates from. That being said, there are tools that can help you here:

salmon (as you mentioned) to quanitfy. Run it with --numGibbsSamples 50 (or higher if your computer has the hardware. Also run with -d for the next tool
terminus is a tool that also originates from the COMBINE-lab, It's relatively experimental, but it applies well to your use case. It will take the salmon output and assign counts to transcript groups. If a set of transcripts can be clearly quantified above a certain confidence, that transcript will stay unique. If it is confounded with other transcripts, those will be placed in a transcript group and their counts will be summarised. More info here

You are right in that TPM, FPKM etc. are not appropriate measures to compare across samples. TMM, or DESeq2's median ratio method (what you get with counts(dds, normalized=TRUE) are a step in the right direction, but are still inappropriate across samples. I'd suggest this often posted blog post in the topic.
There are no perfect solutions here since RNA-Seq is inherently a relative quanitifcation method so this will always be an issue. IMO, use DESeq2::varianceStabilizingTransformation(), DESeq2::rlog() or limma::voom() after terminus and rank your transcripts that way.
Just keep the inherent transcript incertainties in mind.

I did not know about terminus, it looks really interesting. Could you explain why TMM/median ratio method are still inappropriate across samples? I thought they were designed for between sample normalization.

Don't get me wrong, I think both TMM and median-ratio are good enough for your approach, I just don't like to throw around the term normalization lightly. IMO normalization implies aligning dsitributions to allow for statistical testing, which those two don't fully do. If there is a huge disparity in the number and magnitude of expressed genes between conditions, or lots of zero inflation, those methods do not produce a clean normalization either.

If I read you correctly, you want to find the gene with the highest expression?
With RNA-seq data from multiple samples, I'd find this by ranking each gene by TPM/FPKM and taking the mean rank across samples. This method will tend to mitigate any large differences in a gene's expression between samples.
As you are comparing (RNA-seq) expression across genes you must take into account transcript length.


text:
How to incorporate negative controls in DESeq2
We are doing a comparison between two outcomes (positive and negative). We could not have any positive controls as we do not have any "control" data to set as baseline, either from literature or otherwise. No spike-ins have been used. We have some samples which had no resolved outcome, hence we are skeptic of using those as Controls.
We have few negative control samples which contain just the blank culture media with no RNA. We want to remove these signals from the RNA-seq data, but I would like to know the most ideal way to perform this. I would assume just normalizing and subtracting counts wouldn't be the most ideal.
So here's what the coldata looks like -
Sample_ID   Outcome
S1  Positive
S2  Negative
S3  Unknown
S4  Negative
…   …

My previous post has been about the same experiment
Here are the counts for the top 5 genes with high counts in negative samples (RIF1, GOLGB1, C10orf67, MAP3K15, F8):
Geneid  S1  S2  S3  S4  S5  S6  S7  S8  S9  S10 S11 S12 S13 S14 S15 S16 S17 S18 S19 S20 S21
ENSG00000189283.11  2284    2986    4577    2447    1290    1882    1953    1366    2443    3932    2640    1078    491 2919    1248    1627    1805    2094    2229    2789    190
ENSG00000173230.17  1000    1219    1324    1485    722 1142    870 943 1390    1617    1093    564 86  924 778 1059    1144    1248    1169    1828    50
ENSG00000180815.15  641 723 1118    1758    529 865 683 878 1554    2421    1571    414 136 1076    599 924 1098    1199    1405    1327    54
ENSG00000179133.15  369 433 592 488 305 366 335 252 455 2119    1357    118 159 314 200 284 376 333 517 550 17
ENSG00000256637.9   429 491 523 550 250 410 406 359 552 438 339 176 38  357 296 433 517 514 631 779 24

We plan to subtract the median count for these genes from the count matrix of all samples. If the value dips to negative, they will be replaced by 0. Is there a better way to handle background from negative controls? After all these signals are present in all samples, so I'm inclined to believe these signals are noise from the media itself, rather than anything biologically significant from the cells.


text:
How to z-transform Fst and -log(p) values for genome wide selection scan?
I saw in the literature a lot genome wide scans showing $F_{ST}$ values in a 'z-transformed' format which seems to be bounded between 0 to infinity.
Examples are found easily searching for 'genome-wide scan $F_{ST}$' and here is a paper where I could find such z-transformed $F_{ST}$:

Identification of genomic variants putatively targeted by selection during dog domestication (see figure 1)
Genomic divergence during feralization reveals both conserved and distinct mechanisms of parallel weediness evolution (see figure 3)

How are the z-transformed $F_{ST}$ values calculated?
The only reference I could find so far is this: Calculate $F_{ST}$to draw Manhattan plot. And this Population analysis of the Korean native duck using whole-genome sequencing data where they just take the ($F_{ST}$ - mean $F_{ST}$)/sd_$F_{ST}$. But this is giving a lot of negative values, so should this be also surrounded by an absolute value?
Then is it possible to use the z-scores to get -log(p-values)?

The standard score can be estimated by using the sample mean and sample standard deviation as estimates of the population values. Just as you wrote it in the question.
The z-score is given by
$ z = \frac{x - \bar{x}}{s}$
where:
$x$ is the sample observation
$\bar{x}$ is the mean of the sample,
$s$ is the standard deviation of the sample.
Negative values mean that particular observation is $z$ standard deviations below the mean value, and negative $z$-scores are very much expected. The process of standardization is a simple linear transformation and thus does not remove or change any degree of skewness in the original data. I guess your data is skew from the begging to values below the mean.
Yes you can compute the exact $p$-values from the $z$-scores and then apply a negative log transform. You can use a $z$ pre-computed table for this. Check this video for an example.

Thanks for this @JRodrigoF this guy (YouTube link) for teaching critical probability is a teaching maestro.

@JRodrigoF is the correct answer. I'll attempt to explain why the z-score is being used, because it is not immediately obvious.
$F_{ST}$ is not trivial for anyone to understand per se. It is extremely useful for assessing migration - but what is the actual index that Ron Fisher devised?
Simulations suggest that $F_{ST}$ = 0.2 is threshold above which which strong population structure is considered. Thus $F_{ST}$ = 0 panmixia and $F_{ST}$ = 1 is zero migration. How to interpret everything in-between is difficult, except for the 0.2 threshold if the theoretical assumptions in the threshold calculation apply to real-world data. Thus the Z-score is an interesting concept.
In other words, the 0.2 threshold simulation suggests $F_{ST}$ is a non-linear index, which is very interesting in understanding it. Non-linear indices are a bit hairy because we all think they are linear, so 0.1 to 0.3 is one increment of migration, 0.3 to 0.5 another increment of migration. This is not true according to the simulation studies.
Z-scores get around this difficult by simply standardising the combined $F_{ST}$ and thats why they are using it.
It becomes particularly interesting if $F_{ST}$ is calculated per allele because then "fitness traits" can then be tracked. The allele doesn't necessarily need to be the actual fitness trait because of molecular "hitchhiking" (old John Maynard Smith terminology). Thus $F_{ST}$ Z-score within a population within a genome is the beginning of a datamining process. So example, the domestication of dogs why it becomes such a big deal, particularly if the comparison is a dog breed that is closest to the wolf genome. The actually traits selected for in dogs can now be identified. It would also operate between genomes, to identify outliers in dog dispersal, but the behavioural (behavioral) modifications that make dogs dogs is very interesting.

The question is are super small probabilities generated to make log meaningful? I'm not sure in this instance probabilities well beyond very highly, highly significant are generated. If they are then a log transformation is useful.
For background the log-scale is used for transforming very small probabilities rather than use e nomenclature. It makes them more understandable for comparative purposes and allows you to do cool things like add and subtract them instead of multiplying and dividing. The smaller the probability value the higher the transformed value ... so in this case the higher the value the more interesting the result.
However, if the field simply does not use log transformation then it's not cool, because no-one will be able to comparatively interpret the transformed values. In phylogenetics it's done all the time, because the probabilities are really small, but anyone can immediately understand the values - because everyone uses this system. A key part of stats is clarity, not necessarily for those outside the field, but certainly for those within the field.

*very minor quibble it was Sewall Wright not Fisher who invented f-statistics

@user438383 okay, but F is for Fisher (the statistic is named after him), so he must have done something with the development of the stat to get his name recognised like this.

Actually these are not Fisher's F statistics (in fact Fisher didn't develop the F-test, it was just named after him), but Sewall Wright's family of F-statistics, or fixation indices. Confusing, I know...

@user438383, ok I dunno. The reality is Fisher invented so much of classical statistics prior Bayes and maximum likelihood, discovering something he didn't invent really surprises me.


text:
“IndexError: list index out of range” in snakemake
I am trying to establish pipeline in snakePipes and while running DNA mapping I got next IndexError error.
(snakePipes) ~ $ DNA-mapping -i /Volumes/HD2/snakePipes_CHIP/ChiPseq -o /Volumes/HD2/snakePipes_CHIP/ChiPseq1 --mapq 5 -j 10 --dedup hg38
Building DAG of jobs...
Using shell: /bin/bash
Provided cluster nodes: 10
Job counts:
    count   jobs
    5   Bowtie2
    5   FASTQ1
    5   FASTQ2

rule origFASTQ2:
    input: /Volumes/HD2/snakePipes_CHIP/ChiPseq/SRR6761499_R2.fastq.gz
    output: originalFASTQ/SRR6761499_R2.fastq.gz
    jobid: 75
    wildcards: sample=SRR6761499

RuleException in line 9 of /Users/v/miniconda3/envs/snakePipes/lib/python3.6/site-packages/snakePipes/shared/rules/FASTQ.snakefile:
IndexError: list index out of range
Error: snakemake returned an error code of 1, so processing is incomplete!

When I am looking in FASTQ.snakefile:
rule origFASTQ1:
      input:
          indir+"/{sample}"+reads[0]+ext
      output:
          "originalFASTQ/{sample}"+reads[0]+".fastq.gz"
      shell:
          "( [ -f {output} ] || ln -s -r {input} {output} )"

line9: rule origFASTQ2:
      input:
          indir+"/{sample}"+reads[1]+ext
      output:
          "originalFASTQ/{sample}"+reads[1]+".fastq.gz"
      shell:
          "( [ -f {output} ] || ln -s -r {input} {output} )"

I do not understand why I see this error. I hope somebody can point me on the source of the error.
When I use --local parameter I do not see this error, however I get another error:
ln: illegal option -- r

usage: ln [-Ffhinsv] source_file [target_file]
       ln [-Ffhinsv] source_file ... target_dir
       link source_file target_file
    output: originalFASTQ/SRR6761495_R1.fastq.gz
    shell:
        ( [ -f originalFASTQ/SRR6761495_R1.fastq.gz ] || ln -s -r /Volumes/HD2/snakePipes_CHIP/ChiPseq/SRR6761495_R1.fastq.gz originalFASTQ/SRR6761495_R1.fastq.gz )
        (one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)

as far as I understand -r option creates symbolic links relative to link location
in /Volumes/HD2/snakePipes_CHIP/ChiPseq I have:
Bowtie2         SRR6761495_R1.fastq.gz  SRR6761495_R2.fastq.gz
command.sh      SRR6761499_R1.fastq.gz  SRR6761499_R2.fastq.gz
ChIPdict.yaml   SRR6761497_R1.fastq.gz  SRR6761497_R2.fastq.gz
deepTools_qc    SRR6761498_R1.fastq.gz  SRR6761498_R2.fastq.gz
fetch.sh        SRR6761502_R1.fastq.gz  SRR6761502_R2.fastq.gz
filtered_bam
Sambamba
rename_files.sh

I am working on macOS Catalina in bash shell and ln -r is valid option (I installed coreutils)
Thanks


text:
Clustering information saved in Seurat object
Seurat V2 had a option to find clustering information saved in object: PrintFindClustersParams(object = pbmc). How can I get the same clustering parameters from objects in Seurat3? Do I need to manually save clustering parameters or they are saved automatically?
Update: the command line object@commands works perfectly for ...$FindClusters. However, if I would like to find information about dimensions used there is an error:
...$FindNeighbors.SCT.pca Error in [[.Seurat(x, i, drop = TRUE) : Cannot find 'FindNeighbors.SCT.pca' in this Seurat object

I am not sure of the output of PrintFindClustersParams, as I haven't use it before, but try with object@commands. It returns a list containing every command and parameters used.
For example
...
$FindNeighbors.SCT.pca
Command: FindNeighbors(seurObj, dims = 1:20, verbose = vrb)
Time: 2019-12-11 12:41:12
reduction : pca
dims : 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
assay : SCT
k.param : 20
compute.SNN : TRUE
prune.SNN : 0.06666667
nn.method : rann
annoy.metric : euclidean
nn.eps : 0
verbose : FALSE
force.recalc : FALSE
do.plot : FALSE
graph.name : SCT_nn SCT_snn

$FindClusters
Command: FindClusters(seurObj, verbose = vrb)
Time: 2019-12-11 12:41:17
graph.name : SCT_snn
modularity.fxn : 1
resolution : 0.8
algorithm : 1
n.start : 10
n.iter : 10
random.seed : 0
group.singletons : TRUE
verbose : FALSE
...

Re: update
I think it depends on how you built the object. FindNeighbors.SCT.pca should be result of running sctransform normalization. If you use the standard way, you may have other names (e.g.FindNeighbors.pca).


text:
Easier / lazy way to convert an existing bash pipeline to nextflow pipeline?
Not going to lie, my head hurts when I read the NextFlow documentation (No dig on the documentation, I think my head is losing its processing capacity lol).
I had written a pipeline in Bash, essentially iterates over a text file containing full paths of R1 and R2 reads, followed by a series of commands and liberal use of gnu parallel to improve efficiency.
Is there any template I can follow, where I can slot in my commands and get things running quickly?
EDIT - The skeleton of pipeline goes as follows (As I'm not sure if I can put our inhouse pipeline in its entirety here) -
Bash code which loops through R1 and R2, aligns with BWA-mem2 and sorts the bam, then using gnu parallel, I essentially follow the steps as outlined in GATK Best Practices (looking to drop BQSR, but the steps remain) to obtain a gvcf file. The steps in Best Practices have been turned to functions, rather than each step being called sequentially.

Please [edit] your question and show us your code. Nobody can help if they don't know what your code actually does. But, that said, every pipeline is different so I really doubt there will ever be a "template".

Hi KarthikNair could you kindly supply the code or else a reproducible example? Otherwise there are a lot of NextFlow questions and answers on Bioinformatics SE, particularly answered by user @Steve, and those would likely provide pointers. To be honest it's not possible to answer the question as it stands in my personal opinion. I've a reasonable understanding of NextFlow, but I'm not an expert.

I have mentioned the brief gist about the pipeline, essentially we are following GATK best practices for variant calling (humans).

If your objective is to convert Bash scripts into Nextflow code, you might find Viash (https://viash.io) a particularly useful open source tool. It automatically generates independent Nextflow modules from each of your scripts, along with the corresponding Docker containers for each module.
Docker containers are generated on the fly based on your config description in Viash.
These modules can be coupled together to build a (modular) Nextflow pipeline.

Thanks for answering the question as written! Super helpful!

Generating Docker containers on the fly is part of the behaviour of NextFlow per se, which is why it's so cool. However, I know zero about Viash. It would be useful if someone could give feedback on Viash.io because it may not possibly work as advertised here .... The question was flagged for hidden issues, i.e. it may not work as advertised, so I would recommend proceeding with care.

The gist would be something like:
fastqs = Channel.of(...some annoying groovy code to make this a list of tuples...)
workflow{
  align(fastqs)
  makeGVCF(align.out.bams)
  ... collect() the output from above joint call them all...
}

process align {
  input:
    tuple val(sampleName), path(reads)
  output:
    tuple val(sampleName), path("${sampleName}.bam"), emit: bams
  """
  bwa ... | samtools sort ...
  """
}

process makeGVCF{
  input:
    tuple val(sampleName), path(bam)
  output:
    path "${sampleName}.gvcf", emit: gvcfs
  """
  ...haplotype caller code...
  """
}

The usual way would be to use either the fromFilePairs factory method, or to read from a samplesheet using the fromPath factory and the splitCsv operator, e.g. Channel.fromPath( params.samplesheet ).splitCsv( header:true ).map { row -> tuple( row.sample, [ row.fastq1, row.fastq2 ] ) }.set { readgroups }. Then align( readgroups ).

I really thank you both for your inputs. But my head started hurting when I reached the end of the comments. I guess I just have to grit my teeth and push through it.

Since you asked about templates, if your task is something rather generic, there is a good chance that a module, or several, that you can use to build your workflow exists in https://nf-co.re/modules
Generally I would recommend taking a look at the source of nf-core pipelines (e.g. https://nf-co.re/sarek for variant calling) to get an idea how to get your samplesheets in.

Nextflow (https://nextflow.io/) is a workflow manager / language used to implement analysis pipelines, nf-core is an open source collection of modules and pipelines written for / in nextflow

Hi @Niklas thank you - I see your point. Its leveraging NextFlow and providing a link to a NextFlow template on GitHub. Got it. FYI I can use NextFlow. Its Docker container on the fly, is neat, but I agree with DevonRyan its Groovy (an obsolete language) I personally think is a bit weird.


text:
Is there any tool for fuzzy sequence matching?
Note: this question was also asked on biostars
are you aware of any tool that is able to perform error-tolerant pattern-matching search on protein FASTA files?
For example, I want to know, which proteins in my fasta file to matchADNG..C.G regexp (which represents ADNGCG pattern). However, I want to be tolerant to matching errors, meaning that I'm good with any protein that differs in any 1 letter in the motif: ADNG..C.D, ADNG..D.G, MDNG..C.G, etc are all good. Running all possible variants through grep is possible, but exponentially long for longer patterns (I usually have 15-20 letters in pattern, and scan ~10e6-7 sequences).
I am aware of the tool agrep (docs. As far as I know, it's not supported anymore, and also does not distinguish letters in which it's allowed to make error. Also, it does not support long enough patterns (with more than 9 errors -- and yes, I tried to recompile it with myself from here and could not get it work). Also, it's not designed for proteins specifically.
I am also aware of ScanProsite and Protein Pattern Find services. The first one is not error-prone, the second one is seems to be a web interface for grep.

As correctly mentioned on biostars, where I cross-posted it, I should add cross reference to that topic: https://www.biostars.org/p/409111/#409142

Thanks for the link. I personally have no problem with cross-posting, especially if you're linking to other sites and acknowledging the source of any gained knowledge that's added here.

Have you tried using blastp with input containing IUPAC ambiguous bases, and post-filter results to match your criteria? I'm uncertain how ambiguous bases are treated for proteins, as opposed to nucleotide searches, but that might be an avenue to explore.

Python solution:
e<=2 will allow up to 2 insertions, 2 deletions, or 2 SNPs in any 2 combinations
import regex
seq_list=["GGGAGGG", "GGGAAGGG", "GGGAAAGGG", "TTTATTT", "TTTAATTT", "TTTAAATTT", "GGGG", "GGG"]
for sequence in seq_list:
    if regex.search("(GGGGGG|TTTTTT){e<=2}", sequence):
        print(sequence)

GGGAGGG
GGGAAGGG
TTTATTT
TTTAATTT
GGGG
Can also include a wildcard
seq_list=["GGGAAAAAGGG", "GGGAAAAATGG", "GGGAAAAATTG", "GGGAAAAATTT"]
for sequence in seq_list:
    if regex.search("(GGG.*GGG|TTTTTT){e<=2}", sequence):
        print(sequence)

GGGAAAAAGGG
GGGAAAAATGG
GGGAAAAATTG

Agreed on the approach, it can be solved with an outer loop forcing s floating wildcard [AGTC] or [A-Z]  across the regex search term

Thanks for the solution.
The problem is, it is unadequately slow -- e.g. matching whole uniprot with that is very slow.
Also, I am not aware of it's computational complexity for python regex fuzzy matching -- are you aware of its description?


text:
Does the `mmcif/pdb` format store atomic coordinates in units of Ångstrom, nanometers or something else?
I was surprised as to how hard it was to find a reference on this among the docs.
Does the Protein Data Bank mmcif/pdb-x archival format use Angstrom to specify cartesian x,y,z coordinates of atoms or something else or abstract units (perhaps defined per each mmcif _data frame?
How is the origin of the grid calculated -- for crystallographers out there?
Thank you

From google I got this as first hit :
PDB_format_1992.pdf :

Hopefully they keept the same Unit of measurement during the data/file format evolution/extension
It seems the case for mmcif Data Item _atom_site.Cartn_x :

Some history here :
The Protein Data Bank: a historical perspective , The Protein Data Bank
AND :
PY  - 1971
DA  - 1971/10/01
TI  - Crystallography: Protein Data Bank
JO  - Nature New Biology
SP  - 223
EP  - 223
VL  - 233
IS  - 42
SN  - 2058-1092
UR  - https://doi.org/10.1038/newbio233223b0
DO  - 10.1038/newbio233223b0
ID  - ref1
ER  -

About the axis origin CHARMM-GUI has a nice introductory lesson to PDB format (see CRYTST1 Record Section) Lesson 1: Introduction to PDB Format
and stolen from the web, free part of an old book :

Thanks, pip. Missed the "angstroms" part when browsing through the spec. The paper and book references are very much appreciated.


text:
How to grep with a list of genes from a gff the GO terms
I have a text file with a list of name of genes (alias names) like this (** added for emphasis):
Solyc00g019980
**Solyc01g005020**
..
.

And i have a gff file like this:
SL4.0ch01   maker_ITAG  gene    46201   59605   .   -   .   ID=gene:Solyc01g005020.3;Alias=**Solyc01g005020**;Name=Solyc01g005020.3;length=13404
SL4.0ch01   maker_ITAG  mRNA    46201   59605   .   -   .   ID=mRNA:Solyc01g005020.3.1;Parent=gene:Solyc01g005020.3;Name=Solyc01g005020.3.1;Note=Salt Overly Sensitive 1;_AED=0.25;_QI=241|1|1|1|0|0|23|311|1151;_eAED=0.25;**Ontology_term=GO:0006812,GO:0015299,GO:0016021,GO:0055085**
SL4.0ch01   maker_ITAG  three_prime_UTR 46201   46511   .   -   .   ID=three_prime_UTR:Solyc01g005020.3.1.0;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    46201   46844   .   -   .   ID=exon:Solyc01g005020.3.1.1;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  CDS 46512   46844   .   -   0   ID=CDS:Solyc01g005020.3.1.1;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    47281   47585   .   -   .   ID=exon:Solyc01g005020.3.1.2;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  CDS 47281   47585   .   -   2   ID=CDS:Solyc01g005020.3.1.2;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    48276   48543   .   -   .   ID=exon:Solyc01g005020.3.1.3;Parent=mRNA:Solyc01g005020.3.1

I would like to grab only the ones i have in my list gene text file (Alias names) and the corresponding "Ontology_term" next to it.
Terminal is more convenient if there is a solution like grep or something.
Thanks in advance
edit:
i would like to get this output, as an example
Solyc01g005070  Ontology_term=GO:0005543
Solyc01g005230  Ontology_term=GO:0008168

i dont mind if this is included "Ontology_term=", i would prefer not to be included though
thanks again
grep -f list.txt data.gff doesn't work for me, because it grabs all the other as well.

Are you trying to grep the Ontology term from the next line in your example? If so, that makes things a little more complicated.

@JimHat your output example does not correspond to your input examples. Solyc01g005070 and Solyc01g005230 are not referred anywhere in your gff example. Thus, it's hard to relate exactly what you want exactly

@Jim Hat
I took a stab at this given the example files you included.
If the 'list.txt' file looks like this:
Solyc00g019980
Solyc01g005020

And the data.gff file looks like this:
SL4.0ch01   maker_ITAG  gene    46201   59605   .   -   .   ID=gene:Solyc01g005020.3;Alias=Solyc01g005020;Name=Solyc01g005020.3;length=13404
SL4.0ch01   maker_ITAG  mRNA    46201   59605   .   -   .   ID=mRNA:Solyc01g005020.3.1;Parent=gene:Solyc01g005020.3;Name=Solyc01g005020.3.1;Note=Salt Overly Sensitive 1;_AED=0.25;_QI=241|1|1|1|0|0|23|311|1151;_eAED=0.25;Ontology_term=GO:0006812,GO:0015299,GO:0016021,GO:0055085
SL4.0ch01   maker_ITAG  three_prime_UTR 46201   46511   .   -   .   ID=three_prime_UTR:Solyc01g005020.3.1.0;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    46201   46844   .   -   .   ID=exon:Solyc01g005020.3.1.1;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  CDS 46512   46844   .   -   0   ID=CDS:Solyc01g005020.3.1.1;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    47281   47585   .   -   .   ID=exon:Solyc01g005020.3.1.2;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  CDS 47281   47585   .   -   2   ID=CDS:Solyc01g005020.3.1.2;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    48276   48543   .   -   .   ID=exon:Solyc01g005020.3.1.3;Parent=mRNA:Solyc01g005020.3.1

(Note I removed the **s, which I assumed were added for emphasis)
You can do the following:
aliasNames=$(grep -f list.txt data.gff | cut -d':' -f2 | cut -d';' -f1)
ontologyTerms=$(grep -f list.txt data.gff | cut -d';' -f 8 | cut -d'=' -f2)
paste <(echo "$aliasNames") <(echo "$ontologyTerms")

The first line filters the file based on your list and then extracts the alias names from the filtered GFF file. The second line filters the GFF and extracts the ontology terms (after the equal sign). And the third line joins the two as columns.
The result that is printed to the screen is:
Solyc01g005020.3
Solyc01g005020.3.1      GO:0006812,GO:0015299,GO:0016021,GO:0055085
Solyc01g005020.3.1.0
Solyc01g005020.3.1.1
Solyc01g005020.3.1.1
Solyc01g005020.3.1.2
Solyc01g005020.3.1.2
Solyc01g005020.3.1.3

Given the examples you included all of the lines of the GFF are included, but not all have gene ontology information.


text:
How to introduce point mutations in a protein structure?
I want to introduce point mutations to a protein PDB file and automatically predict the 3D coordinates of a new amino acid residue in the mutated protein.
I used interactive Schrodinger's Maestro tool to mutate residues and minimize them to avoid clashes, however this time i want to automate the task and use some open source tools.
Is there any free software for this task?

PDBFixer and pyMutate are simple, lightweight, free and open-source. Both don't do any optimizing, though, so you may have to run some molecular dynamics afterward to get what would be the right "orientation".
VMD would work for sure, I believe PyMOL too. They are both graphical interfaces but you can start them in text only mode (or just import pymol in a python script and use the Python API). They are free for particulars and academics, but given their renown you may already be using them.
If you're looking for something a tad more complete and that may adapt the side-chain to the surrounding environment I'd go with Modeller, which provides two example scripts to do point mutations (here and there).
If you don't want to just get the side-chain of this amino-acid right but actually strive to see some general backbone adaptation too then use PyRosetta. It can probably do what all the others do, too, but given it is 1 GB big it may be overkill. Any of the package coupled with molecular dynamics would also do the trick.

Dear @O.Laprevote, thank you for your valuable answer! I actually tried running Rosetta's fixbb tool to mutate the truncated pdb residue with all missing atoms but one, however apparently it just ignores such residues, while mutating the truncated residues was my very task.

Modeller did the job very well for me. I am accepting this answer!


text:
hisat2 --rna-strandness option and downstream htseq-count analysis
I've got some doubts on the hisat2 --rna-strandness option and its output for downstream analysis. Is it expected to see a difference in alignment and counting results given the default usage of --rna-strandness in hisat2 followed by htseq-count -s reverse on a strand-specific assay?
Please see below.
I understand that the --rna-strandness option produces an XS tag to indicate where a transcript is from (on the + or - strand) for downstream transcriptome assembly analysis. I have a paired-end stranded sequencing library that was aligned to the genome using hisat2 without specifying the --rna-strandness (in other words, the default unstranded was the usage). Following this, the reads were assigned to genes using htseq-count and this time -s reverse was specified given the strand-specific sequencing assay type.
Since --rna-strandness is for transcriptome assembly using the XS tags generated and htseq does not use XS tags for counting, I presume there should be no practical impact from the above. Could you also shed light on this? in case I may have been overlooking other facts of the usages of the tools.
To help verify the above, I re-aligned and counted the reads from 2 samples by switching on --rna-strandness RF in hisat2. I attach the alignment and count features info. below for assessment.
Overall alignment rate of Sample 1:
94.52% (--rna-strandness RF)
94.12% (--rna-strandness unstranded)

Overall alignment rate of Sample 2:
94.57% (--rna-strandness RF)
94.15% (--rna-strandness unstranded)

Feature counts of Sample 1 (following --rna-strandness RF + -s reverse):
__no_feature        6327294
__ambiguous     2954776
__too_low_aQual     3784481
__not_aligned       688856
__alignment_not_unique      4858182

Feature counts of Sample 1 (following --rna-strandness unstranded + -s reverse):
__no_feature        6291151
__ambiguous     2911298
__too_low_aQual     4075017
__not_aligned       754400
__alignment_not_unique      16136045

Feature counts of Sample 2 (following --rna-strandness RF + -s reverse):
__no_feature        5417882
__ambiguous     1708510
__too_low_aQual     3532352
__not_aligned       564596
__alignment_not_unique      2859501

Feature counts of Sample 2 (following --rna-strandness unstranded + -s reverse):
__no_feature        5359434
__ambiguous     1676091
__too_low_aQual     3813344
__not_aligned       623122
__alignment_not_unique      2891792

These results look comparable to me across pipelines.
Thanks
Guan

Maybe this post can help https://bioinformatics.stackexchange.com/questions/4074/hisat2-which-option-should-mention-for-strand-specific-library-read/4075#4075

If you reran the command with the correct settings, just leave it at that.  (It is not at all clear to me that strandedness rf is correct)
If you want people to tell you if you ran the commands right, you need to put down what commands you used.

Sample 2 looks comparable but Sample 1 shows a drastic reduction in unique alignments.
Your strandedness param is going to dictate how your alignments are done in the sense that discordant pairs will either be detected or not.
I would identify the strandedness of your library and use the appropriate respective param.


text:
Are there any Sulfur Dioxide molecules in this structure?
Does the protein structure shown in the JSmol link have any Sulfur Dioxide molecules?

In the picture, the Sulfur element is yellow and Oxygen is red (see here)
https://www.viprbrc.org/brc/structureOperation.spg?accession=1P9S&decorator=corona&context=1583712525175
Link to mmcif file (I have the VS Code plugin installed):
https://files.rcsb.org/view/1P9S.cif

The answer, as provided in the comments by gilleain is No. The commenter wanted to know why I asked this. The reason is I heard sulfur dioxide as silent speech when thinking about antibodies that’s why I asked.
I searched the file https://files.rcsb.org/view/1P9S.cif for the SO2 molecule and didn't find it. I got confirmation in the comments of the question.
Would anyone be interested in an app that takes a .cif file of a molecule and allows you to search the file for any type of element or molecule in it? Would that help anyone with their research?

Sorry, I do not see the advantage. All three PDB databases allow ligand search by name, residue name or SMILES in the browser. The PDBe has a decent API (https://www.ebi.ac.uk/pdbe/pdbe-rest-api), which gives info of which ligands are present, without downloading the structure. Furthermore, most instances of elements beyond the 6 of life are crystallisation artefacts, so it is way safer to search for enzymes that are expected to have more exotic elements, e.g. tungsten in a W-enzyme as seen in PDB:2W7Z.


text:
How many single copy gene are needed to construct the phylogenetic tree?
I ran the orthomcl program with genomes of the whole streptomyces genus(from ncbi) as input, but I only got 22 orthologous groups(streptomyces's genome size vary from 2M to 15M, which shows the great difference among them). And the analyst from the bioinformatics company told me that my result wasn't reasonable if I try to construct the phylogenetic tree with the 22 orthologous groups I got. He also told me, I'd better use almost 1000 orthologous groups to build the tree. My question is, how to know my orthologous groups are enough to construct a tree?
The program was run during my summer vacation(a long time ago...), but I still can't realize why "1000 orthologous groups are needed". And the analyst form company just tell me that it was the experience number. But when I read some article about comparative genomics, I see some papers are used just about 10 single-copy genes to build the tree, why???

Thanks, your research objectives are complicated and difficult to answer in a question and answer format, even for a phylogeneticist with experience in bacteria. The other issue is language, because potential typing errors can change the meaning of what you are doing.  Lets keep trying, I'm sure we'll get there.

Okay I'm not really sure so I've bountied your question. Please supply as much information as possible  about your 22 contigs, lets see if someone else could anwer this.

The one thing that strikes me is that OrthoMCL was written under David Roos supervision, which really means it was written with Toxoplasma gondii in mind, i.e. a eukaryotic parasite with notable levels of gene duplication. They tested it against Plasmodium falciparum, but that isn't surprising because both are members of the Apicomplexa. Personally I would suggest looking at a tool more orientated towards bacteria, albeit I'd like to see the output of OrthoMCL.
Prokka is a tool that might give better results, here, which is designed for bacteria, but personally I would extract the orthologes/homologues directly if the primary purpose is to contruct a phylogeny, i.e. I would parse annotated genes which are not "hypothetical".

Just wanting to add that I've had good experiences with using Orthofinder on bacteria, and it's more user-friendly than I'd have thought possible for something that uses a command line (I'm leaving this as a comment because I've never used OrthoMCL).

I would say, 120 markers is enough for phylogenetic analysis.
see https://gtdb.ecogenomic.org/


text:
Troubleshooting E. coli Strain Identification: Identifying Transposon Insertion Sites with Paired-End Illumina Data
I have a problem that needs your help.
I am working on a bacterial genomics project working with E. coli. The lab did illumina sequencing of many strains but I have some problems identifying the WT version (MG1655) from the recA, and ruvAB recombination gene mutants. The mutations were done by  gene interruption 20+ years ago (I suppose using tn5). While they genetically verified that the genes do not work (also many years ago), they never checked what is the sequence of those mutations, for example, making a PCR and Sanger sequencing it. Which lead to my problem.
I am analyzing the sequencing results of the WT, recA, and ruvAB strains and the results makes more sense if the strains were in a different order: recA, ruvAB, and WT. Therefore, I need to verify that the strains are in a different order. Sequencing was pair-ended and I got about 20 million 50bp reads. Mapping the reads against the reference genome does not show a decrease on the depth of sequencing on any region either of the recA and ruvAB (because the gene got interrupted by a transposon) or at least not enough to consider that the transposon might be inserted on that gene. Which leaves me wondering how can I use the Illumina sequencing reads to verify which mutant does that genome has.
Is it possible to map against the reference genome and ask not to get the reads that map into it, but the ones that don't? That way, the ones I am left with should include the transposon. Then, maybe, assemble them de-novo. Check that the transposon is a Tn5 and what other nucleotides are left at the ends of the tn5 sequence. Finally, check those ends against either the recA or ruvAB gene.
Does this approach makes sense and can be done? If so, how? Does someone has a better idea on how to solve this problem?
Many thanks!

If you're asking about negative selection from mapped reads, yes, that can be done. There are many ways to do it, but the easiest is to use a mapper to create a BAM file containing all reads, then apply a filter on that file to extract the unmapped reads.
An example process would be as follows:
bowtie2-build reference.fasta reference.fasta
bowtie2 -x reference.fasta -U reads.fq.gz | samtools sort > reads_vs_ref.bam
samtools view -b -f 0x04 reads_vs_ref.bam > unmapped_reads_vs_ref.bam


text:
Filtering genes from cuffdiff results
This question has also been asked on Biostars
I have run cuffdiff (with statistics turned ON) to compare two groups of samples: Control group and Late AD group.
This is the command I ran to be precise:
cuffdiff -L Control,AD_Late_Braak -p 8 --total-hits-norm --frag-bias-correct ../References/ensembl.GRCh38.99.fa --multi-read-correct --library-norm-method quartile ../References/Homo_sapiens.GRCh38.99.chr.gtf Early_Braak_Control_1,Early_Braak_Control_2,Early_Braak_Control_3 Late_Braak_Sample_1,Late_Braak_Sample_2,Late_Braak_Sample_3

I'm looking at the output from cuffdiff and I see a gene_exp.diff file which contains the results of the differential expression testing. I want to know what is the best way to filter the results of this gene_exp.diff file so as to restrict the number of genes that are up-regulated and down-regulated to a list between ~50-850.
P.S My thoughts are to tweak the p-values and log2 fold change values, but seems like a trail-and-error method, so I was wondering if there was a more "formal" method/approach?

To filter your gene_exp.diff either you export it to excel and then filter because it is columns separated by tabs or use the command (awk and head -n) to filter and have the number of lines you want if you master Linux


text:
How to create this function for lapply?
I have two lists of dataframes. List x and list my_list. Each having four dataframes and their names are same for the two lists, i.e., bray-curtis, Chao1, FB, and Shannon.
The dataframes within the two lists have same element names, but they contain different information. This is because list x has been generated after some operations on list my_list.
Now I want to write a function and run that on list x with `lapply. the command I want to write as a function is
forest(x$Chao1, slab = paste(my_list$Chao1$StudyID))

Here, the first Chao1 is coming from list x while the second  Chao1 is coming from the list my_list. It works for a single dataframe but I am unable to write the codes with lapply so that it runs across all the dataframes of the lists.
The desired outputs are four plots from four dataframes of list x.
Can anyone please help me?
Update: @gringer's approach works, but it saves all the plots in a single list. When I tried to save the images with following command it just saved the last image:
tiff(filename = "~/Desktop/image.tiff", width = 30, height = 30, units = "cm", res = 1200)

Is there a way to save all them individually?
Thanks

I am still quite confused. You mentioned that the names in my_list and x are the same, why do you need the names from my_list then? Also note you can [edit] your question and add all these details. Try for example, showing bits of the data.frames in the two lists, or something that would make us understand why you need two lists, not one.

Assuming the lists are named, and you want to fetch the element with the same name from each list, I'm guessing something like this would work (or at least get you closer to a solution):
v <- lapply(names(x), function(var){
   forest(x[[var]], slab=(my_list[[var]])$StudyID)
})

When creating images that don't have a multi-page format (i.e. anything except PDF), and not changing the image name with each save, you need to add in a format specifier to tell R what to do with new created images. For example, this will add a three digit zero-padded number for each new TIFF image:
tiff(filename = "~/Desktop/image_%03d.tiff", width = 30, height = 30, units = "cm", res = 1200)

Thanks gringer, it really worked. Being new in R, I would like to ask how does this line of code work? I mean, why I have to put double brackets [[]] there? Thanks


text:
How to import data from cell ranger to R (Seurat)?
I will have some scRNA-seq data. The goal of the experiment will be to see if there is any difference in gene expression between treatment groups using the package Seurat from R.
I have read a tutorial how to do the analyze, but this tutorial does not explain how to import data.
I am a total beginner with R, I searched the internet,but the majority of the explanations related to csv/excel/txt... common files.
Data are from Cell ranger and spread in 3 files with following file extensions : .tsv and .mtx (barcodes.tsv, genes.tsv and matrix.mtx).
Collaborators ran Cell Ranger and gave these cell ranger output files : barcodes.tsv, genes.tsv and matrix.mtx.
Can someone give me the code to import these kind of data to R ?
Thank you in advance.
update: By "directory" you means "folder" on my computer for example ? So I need to save the 3 files under one folder and then to put name of the folder into brackets after the Read10x function ? I am a total beginner with R (I have used to R commander for statistical tests but that's it...).
update 2: I successed to import data into Seurat with the explanation of StupidWolf. I do not know C language ( I do not know almost any language ..). My system crashed a little further in the analyze after importing, not enough RAM... so now I am learning Linux commands to switch on the cluster.

Answer from StupidWolf, converted from comment:
There is a function call Read10X in seurat, you can see how to use it with the example in https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/rkit
You place the three files into a directory / folder, then you specify Read10X()... hey how big is your data? you might need to run this on a machine with enough RAM.
You can try downloading the data following instructions here:
https://satijalab.org/seurat/v3.1/pbmc3k_tutorial.html

Have a look here and the entire tutorial:
https://github.com/hbctraining/scRNA-seq/blob/master/lessons/03_SC_quality_control-setup.md

To convert data from one format into another, you may find some software tool sometimes. But mostly no such tool can be found to satisfy all the requirements that we want.
In such situations, I recommend that you write a C program that reads data from the input file and writes it into the output file. In this way, you have complete flexibility so that you can determine the exact output format and what part of data to ignore or include. To do that, you should be familiar with some concepts in C language such as arrays, loops, strings, and file operations.
For an example of such conversions, you can get ideas from the following page:
https://stackoverflow.com/questions/55347322/c-program-to-convert-a-text-file-into-a-csv-file
Of course, such a converter program can be written in other languages too such as C# and Python.
The program can be developed in R language too. You need to be familiar with import functions, loops, strings, and arrays in this language.
We cannot expect that custom tools exist for all our file convertions. Every bioinformatician must be skilled in a standard programming language such as C, C#, Java, Python, or R. Then, he/she will be able to develop any custom program that reads an input file and generates a custom output file.

It's a rather onerous workaround to write C programs as interfaces between R packages and data (R generally has excellent data import functionality).


text:
Theoretical limit of human genome compression
How small can a compressed file containing the human genome be?
I'm aware that this question cannot actually be answered, since it is asking for the Kolmogorov complexity of the human genome, which is not computable. So reformulating:
What is the largest known lower bound that we can put on the size of a file containing the human genome?
An example answer could be: on average we know that the human genome contain at least this amount X of randomness, therefore it is not possible to do better than Y bytes.
I'm agnostic on which human genome, but I would say without polymorphisms to make it simpler.
The compressed sequence would need to be stand alone, without using a reference (I asked a question on that topic here: Compressing the human genome to few megabytes).
There should be a program able to compress and decompress the genome in a lossless way. Apart from that, any kind of modelling and coding is allowed.

For the most modern CHM13 assembly (v2.0; Complete T2T reconstruction of a human genome; https://github.com/marbl/CHM13), the best known compression without any use of a reference is 567 MB (1.45 bits/bp; 1.80 bits/distinct canonical 31-mer).
Upper estimates of genome compressibility can be made using standard general or specialized compressors (see, e.g., this comparison). In my experience, xz can achieve great compression ratios on genomic sequences, especially when used with the maximum compression level (-9) and a single thread (-T1). The compression performance can be further boosted by using a large dictionary (-e --lzma2=preset=9,dict=512MiB). Importantly, the original FASTA file needs to be cleaned – linebreaks hamper compression as they act as noise; the reformatting can be done by seqtk seq -U (which also converts all nucleotides to the upper case). It the case of the CHM13v2 assembly, xz compresses the human genome to 634 MB.
However, even better performance can be achieved with specialized compressors, especially if they have particular support for reverse complements and mutational models. Among the best ones is GeCo3, which has been shown to compress CHM13v2, with a particular combination of parameters, to 567 MB 3.
Lower estimates are impossible in general (we can't reject the existence of a program that would create the whole sequence from small amount of starting conditions), but can be done for particular compression techniques, e.g., particular variants of dictionary compressors. However, this leads to intractable optimization problems (see e.g. the problem of minimal string attractors).
A useful auxiliary measure for estimating the compressibility of single genomes is the number of distinct k-mers for a "reasonably" selected k (e.g., k=31). The lower bound on the compressibility of k-mer sets can then be roughly approximated by 2 bits per distinct k-mer, however, this is limited by the embedded structural assumptions, e.g., that single genomes are in principle random strings with occasional (possibly distant) repeats. For instance, the CHM13v2 human genome assembly contains 2,512,390,070 distinct canonical 31-mers, which would then correspond to 628 MB assuming 2bits/kmers, while the best available specialized compressors can compress it down to 567 MB, i.e., by approx 10% better.
Updates

May 27, 2022 – updated the upper bound based on the latest results with CHM13v2 (xz & GeCo3), changed the lower bound wording to an auxiliary estimate as the assumptions have shown not to hold

There's a benchmark that compared dozens of compressors on DNA data, including human genome: http://kirr.dyndns.org/sequence-compression-benchmark/
In particular, this table shows file size of compressed human genome, produced by the strongest setting of each compressor: http://kirr.dyndns.org/sequence-compression-benchmark/?d=Homo+sapiens+GCA_000001405.28+%283.31+GB%29&cs=1&cg=1&com=yes&src=all&nt=4&only-best=1&bn=1&bm=size&sm=same&tm0=name&tm1=size&button=Show+table
The best result is currently 657 MB on that particular genome build (GCA_000001405.28), when compressed with fastqz.
By the way, xz (mentioned in another answer) is currently number 20 in that table. (Though the genome is un-modified straight from the database, not the "cleaned" form.)
Disclosure: I worked on that benchmark.
Note: What is being compressed here is a FASTA file containing a particular assembly of a human genome. Compared to an actual genome it has this extra content: Sequence names and brakes between contigs, end of line characters, N characters representing unknown bases, and soft-masked repeats (represented as lowercase). On the other hand, it also misses some information from the actual genome: 1. It's an unphased haploid sequence, while the real genome is diploid (two copies of each chromosome). 2. It includes no epigenetic information (DNA methylation, histone modifications), nucleosome positioning, 3D structure of packed DNA, etc. So, it's just a digital model of a genome with a certain degree of accuracy, that is useful,  simple, and relatively affordable to obtain. It's a format commonly used for representing genomes in databases.

I quickly tested xz on this build (GCA_000001405.28) in the same way as in my reply (seqtk seq -U; xz -9 -T1) and got 667 MB, which is very close to fastqz-slow in the referenced table (657 MB). With general compressors, it's important to remove eols in the sequences, as well as to convert all bases to the upper case. Note that none of these steps modify the genome in any way – fastas can always be converted to a fixed width, and the upper/lower case info are usually just masked repeats.

As the table above says 770 MB (i.e., +103 MB) for xz with similar parameters (-e9), I assume the file probably wasn't in a form that would facilitate compression using general compressors, and the observed gain of specialized compressors came mainly from removing the eols/upper case conversion.

An additional experiment: I verified if xz could outcompete even fastqz-slow, and it indeed can. GCA_000001405.28 can be compressed down to 652 MB by 'xz -T1 -e --lzma2=preset=9,dict=512MiB' (again, after 'seqtk seq -U').

@Kamil S Jaron, cleaning is one reason, another is that it's just different genome build. Though as you see Karel obtain similar number with his method.

@Karel Břinda, by the way, it's good idea to tweak xz's dict (dict=512MiB), I should benchmark this setting too. Thanks.

Thanks both, although you are right Kirr, that in real-life situation you don't want to lose information about softmasked bases, the OP asked about the number of bites to store the human genome, therefore Karel's approach is completely valid (human genome does not have really softmasked bases, that's just our way how to add there one more layer of information).
Both of you, would you like to update the answers? I think this discussion is really interesting.

@Kamil S Jaron, I updated the answer. You're right that it's just one way of representing a genome, and also many genomes don't include soft-mask (e.g. prokaryote genomes).