# Question
Title: 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

```

# Answer
> 1 votes
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 :)

# Answer
> 1 votes
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()

```

# Answer
> 1 votes
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))

```

---
Tags: r, ggplot2
---


# Question
Title: How to represent pathways as mathematical graphs?
I would like to construct a mathematical graph: $G = (E,V, w(E))$ a graph with vertices, edges with weight on them based on a KEGG pathway.

Is there a way to do that with existing tools? If not I will start to implement a Python module for that. My main question is what is the mathematical graph representation of pathways? Assuming the possible node types and the different connections?

I would like to transform a pathway to a graph because then one can run min-cut-max-flow algorithms and stuff like that as part of my pathway analysis.

# Answer
Pathways are generally directed graphs, which can be implemented in python using the `networkx` package. Note that there will be cycles in graphs.

Please note that there are large numbers of pathways analysis packages that already exist (both python and R packages but also just webservers), so I strongly recommend searching to see if something already exists before wasting time and implementing it yourself.
> 1 votes

---
Tags: python, pathway, kegg
---


# Question
Title: Ancestry of the coronavirus 2019-nCov, WuHan city, China
In one of the answers to another question about the corona virus a link was given to this phylogenetic analysis of the virus.

Can somebody assist a non-bio type here? It seems to show that the current corona virus split from a virus in bats. And that the same ancestral virus is also ancestral to several other viruses, some in humans and some in bats.

Is my understanding basically correct?

# Answer
> 5 votes
**Overview** The central focus of the tree is to highlight the key biological concern of the new coronavirus, 2019-nCov. The key concern is the genetic similarities to SARS epidemic, and relates to the SARS receptor.

**SARS background** SARS is endemic in bats (your BioRxiv tree partly shows that and this tree definitely shows it) and in the 2002 epidemic infected civet cats which then infected humans. More importantly it is then transmitted from one human to another. This is a consequence of the SARS receptor been able to exploit the cellular receptors in the respiratory tract of bats, civet cats and humans and use this tissue as a site of replication. The concern is exacerbated because the tree here shows SARS independently infected humans on two separate occasions, suggesting the cross-over has an underlying genetic basis, which is part of the "SARS lineage".

**BioRxiv 2019-nCov tree** Your tree shows the 2019-nCov has a recent common ancestor with SARS in comparison to the rest of the betacoronaviruses. This therefore provides some circumstantial evidence the receptor mechanism and ability to frequently crossover from bats and ultimately transmit human to human could be shared with SARS. At an evolutionary level (which is what ANY tree is supposed to reconstruct) it raises the question which common ancestor of SARS did the "SARS receptor" originate. I'd need to draw a diagram to better demonstrate this point, but I hope you get the idea.

**The "influenza receptor"** The analogy is with influenza virus and the entire epidemiology of influenza hinges around is sialic acid receptor of influenza and its ability to bind on to the cellular receptors in the upper/lower respiratory tract of birds, pigs and humans. The hypothesised mechanism is called the mixing vessel theory and is the classic epidemiological understanding of how new pathogenic influenza pandemics occur. If you replace "birds and pigs" with e.g. "bats, civet cats and humans" you get the idea why this could be scary. We don't the intemediate host of 2019-nCov ... but I speculate there must be one, unless eating Cov-infected bats is common in China.

**Technical details of the BioRxiv tree** The tree is a nicely diverse selection of the beta-coronaviruses. The authors have rooted the tree using the outgroups of delacoronaviruses and gammacoronaviruses, so it is a good robust selection of outgroups which can be used to correctly identify the direct of evolution of betacoronavirus divergence. In tree theory (its formal name is phylogenetics theory), extensive rooting is good and minimises artefacts.

One of the earliest members of the betacoronaviruses to diverge is MERS (Middle Eastern respiratory syndrome), which form a single "clade" (all viruses share a unique common ancestor) and these represent around 50% amino acid divergence from SARS/2019-nCov. The selection of MERS in the tree from both camels and humans looks good. The selection of other betacoronaviruses looks great, I wasn't aware of the "ruminant clade" at all involving buffalo, cow etc .. infections and there was an associated human infection. There are loads of bat isolates of the betacoronaviruses throughout the tree, but we eventually arrive at the "SARS clade". The authors show that 2019-nCov is an outgroup to the SARS clade and shows a close relationship to one (BioRxiv tree) or two (this tree) bat isolates. Looking at the precise SARS clade (better shown in this tree) we can see loads of bat virus associating with SARS lineages. We therefore assume the reservoir to SARS and likely 2019-nCov is bat, moreover that the single ancestor to both viruses was a bat (it is called a parsimonious hypothesis).

The one thing the BioRXiv tree omits, for example with regards this tree, is the diversity of SARS and in particularly the two independent origins of SARS, which is a weakness of their analysis, particular if this fed into downstream analysis. It is not to say the authors were wrong, but it was uncool.

**Word of caution** my understanding is the divergence between SARS and 2010-nCov is around 15% and this is a quite a large amount of genetic divergence, even if both viruses share a most recent common-ancestor, form part of the same clade and receptor peptide motifs. Nevertheless there is sufficient amino acid divergence to generate notable differences in epidemiology, clinical symptom and transmission.

**MERS, SARS and 2019-nCov mortality rates** It is worth noting that within the betacoronaviruses the mortality rates between difference clinically important viruses is very different. MERS has a mortality rate of 40-60%, SARS is around 10% but 2019-nCov is \<2.5%. The 2019-nCov mortality rate is still important given it has infected more and could potentially infect alot more people than SARS.

---
Tags: sequence-alignment, fasta, phylogenetics, phylogeny, covid-19
---


title: Converting values to string

text:
# Question
Title: Converting values to string
I have a data frame

```
> head(custom.cn.data)
      Gene       Sample_name         CN
1: DDX11L1          F02         0.6788844
2: DDX11L1          C02        -0.9924524
3: DDX11L1          A07        -2.3833559
>

```

I want to convert values in CN column \< 0 to Del and \>0 to Amp

But I am getting error

```
> a %>% mutate(q = replace(CN, which(CN<0), Del))
Error in replace(CN, which(CN < 0), Del) : object 'Del' not found

> a$CN <- replace(a$CN, which(a$CN < 0), Del)
Error in replace(a$CN, which(a$CN < 0), Del) : object 'Del' not found

```

Any help?

# Answer
```
custom.cn.data$new_column <- ifelse(custom.cn.data$CN < 0, "Del", "Amp")

```

In this scenario `CN == 0` would be `Amp`
> 6 votes

# Answer
The following code would create a new column with the desired values, 0 values will have `NA` in this new column.

```
custom.cn.data[custom.cn.data$CN < 0, "CN_string"] <- "Del"
custom.cn.data[custom.cn.data$CN > 0, "CN_string"] <- "Amp"

```
> 4 votes

# Answer
Try to stay within one package (paradigm) for one problem. This will make the code more readable, easier to debug.

At the moment you have data as a *data.table* object from **data.table** package, then trying to add a new column using *mutate* from **dplyr** package, and within *mutate* using **base** *replace* function.

Solutions:

**data.table**

```
# update by row index
custom.cn.data[ CN < 0, NewColumn := "Del" ]
custom.cn.data[ CN >= 0, NewColumn := "Amp" ]

```

or use *ifelse*

```
custom.cn.data[, NewColumn := ifelse(CN < 0, "Del", "Amp") ]

```

**base**

see @DavyCat's and @haci's answers.

**dplyr**

Note, this will change the *data.table* into *data.frame*.

```
# ifelse
custom.cn.data <- custom.cn.data %>%
  mutate(NewColumn = ifelse(CN < 0, "Del", "Amp"))

```

or use *case\_when*

```
custom.cn.data <- custom.cn.data %>%
  mutate(NewColumn = case_when(
    CN < 0 ~ "Del",
    CN >= 0 ~ "Amp"))

```
> 2 votes

---
Tags: r
---


# Question
Title: samtools / bamUtil | Meaning of <wildcard> as Reference Name
I have been working with several tools on bam files recently and I am not sure how I should interpret some of the outputs.

```
samtools idxstats mapped.sorted.bam

```

```
Chr1    15096   0   0
Chr2    33397   0   0
Chr3    43888   41  0
Chr4    20819   1   0
*   0   0   34

```

Documentation says third and fourth column give "mapped reads" and "unmapped reads". So I interpreted the last line (`*`) as the number of unmapped reads as this is the only line with a non 0 value at this field.

Then I used

```
bam splitChromosome --in mapped.sorted.bam --out test/aln.

```

```
Reference Name: chr6 has 6 records
Reference Name: Chr7 has 1 records
Reference Name: Chr3 has 41 records
Reference Name: Chr4 has 1 records
Reference Name: * has 34 records

```

And looking at the bams outputted (one per reference sequence) in `test/` I find an `aln.UnknownChrom.bam` that I am not sure where it comes from. I actually found 97 unique reads ID in it, where I was waiting for 34 or less.

So my **first question** is : Does "\*" actually means "all reads that didn't mapp" and this information is stored in way I didn't know about in BAM files ?

And my **second question** would be : Anybody knows where this `aln.UnknownChrom.bam` comes from when using `splitChromosome` ? It is a source of error in a pipeline I am building and I'd like to get rid of it, if these are not actual mapping informations.

# Answer
> 2 votes
As per the SAM/BAM format specifications:

> "RNAME: Reference sequence NAME of the alignment. If @SQheader lines are present, RNAME(if not ‘\*’) must be present in one of the SQ-SN tag. **An unmapped segment without coordinate has a ‘\*’ at this field.** However, an unmapped segment may also have an ordinary coordinate such that it can be placed at a desired position after sorting. If RNAME is ‘\*’, no assumptions can be made about POS and CIGAR." -- http://samtools.github.io/hts-specs/SAMv1.pdf, emphasis mine

So, yes, these are the unmapped reads.

I haven't used this `splitChromosome` utility before myself, but as it splits the BAM file into one BAM file per chromosome, I suspect it simply also creates a BAM file for the unknown (`*`) chromosome. Why the numbers don't match I can't tell.

How to deal with that file in your pipeline will depend on what the pipeline is meant to do. But if you aren't interested in these reads, I suppose you could simply remove it before continuing to the next step.

---
Tags: sam, samtools
---


# Question
Title: Old versions of the reference genome and dbSNP
For benchmark purposes I’m trying to find an old version of the human reference genome (pre-GRCh37) and a matching version of dbSNP.

Unfortunately it appears as though older versions of dbSNP aren’t archived, or at least these archives are not publicly accessible.

In particular, the NCBI FTP archive is fairly sparse. The only pre-GRCh37 version it contains is ostensibly dbSNP build 124 from January 2005, for hg17 (NCBI build 35). Unfortunately the corresponding FTP folder is actually *empty*. dbSNP build 138 has been “removed to recover storage”, whatever that means, and I’m unclear which version of the human reference genome it refers to (its official page does not contain this rather crucial bit of information).

Is there another authoritative source of this data?

# Answer
> 2 votes
The UCSC Genome Browser downloads page has old human genome assembly versions way back to the first human genome ever published (by UCSC in 2000). http://hgdownload.soe.ucsc.edu/downloads.html#human

We also have a variety of archived dbSNP data up to dbSNP141 on hg38, dbSNP 131 on hg37, dbSNP126 on hg18: http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/snp141.txt.gz http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/snp138.txt.gz http://hgdownload.soe.ucsc.edu/goldenPath/hg18/database/snp126.txt.gz

They're in a custom table format, which can be turned into VCF manually.

You can also get all of this data from Table Browser and add filters, etc.: http://genome.ucsc.edu/cgi-bin/hgTables

---
Tags: database, public-databases, variant-calling, reference-genome
---


# Question
Title: How to add an increment to a repeated strings
This is my data dput(METADATA)

```
structure(list(Run = c("SRR2920466", "SRR2920477", "SRR2920478",
"SRR2920506", "SRR2920507", "SRR2920532", "SRR2920531", "SRR2920469",
"SRR2920470", "SRR2920481", "SRR2920482", "SRR2920500", "SRR2920501",
"SRR2920536", "SRR2920537", "SRR2920471", "SRR2920472", "SRR2920483",
"SRR2920484", "SRR2920505", "SRR2920538", "SRR2920539", "SRR2920475",
"SRR2920476", "SRR2920487", "SRR2920488", "SRR2920542", "SRR2920543",
"SRR2920555", "SRR2920559", "SRR2920570", "SRR2920573", "SRR2920578",
"SRR2920583", "SRR2920586", "SRR2920594", "SRR2920546", "SRR2920547",
"SRR2920548", "SRR2920549", "SRR2920550", "SRR2920551", "SRR2920553",
"SRR2920558", "SRR2920563", "SRR2920565", "SRR2920567", "SRR2920569",
"SRR2920572", "SRR2920575", "SRR2920577", "SRR2920580", "SRR2920582",
"SRR2920585", "SRR2920589", "SRR2920590", "SRR2920591", "SRR2920593",
"SRR2920554", "SRR2920552", "SRR2920556", "SRR2920557", "SRR2920560",
"SRR2920561", "SRR2920562", "SRR2920564", "SRR2920566", "SRR2920568",
"SRR2920571", "SRR2920574", "SRR2920576", "SRR2920581", "SRR2920584",
"SRR2920587", "SRR2920588", "SRR2920592", "SRR2920595", "SRR2920579"
), source_name = c("HSC", "HSC", "HSC", "HSC", "HSC", "HSC",
"HSC", "CMP", "CMP", "CMP", "CMP", "CMP", "CMP", "CMP", "CMP",
"GMP", "GMP", "GMP", "GMP", "GMP", "GMP", "GMP", "Mono", "Mono",
"Mono", "Mono", "Mono", "Mono", "LSC", "LSC", "LSC", "LSC", "LSC",
"LSC", "LSC", "LSC", "Blast", "Blast", "Blast", "Blast", "Blast",
"Blast", "Blast", "Blast", "Blast", "Blast", "Blast", "Blast",
"Blast", "Blast", "Blast", "Blast", "Blast", "Blast", "Blast",
"Blast", "Blast", "Blast", "Blast", "pHSC", "pHSC", "pHSC", "pHSC",
"pHSC", "pHSC", "pHSC", "pHSC", "pHSC", "pHSC", "pHSC", "pHSC",
"pHSC", "pHSC", "pHSC", "pHSC", "pHSC", "pHSC", "pHSC")), class = c("spec_tbl_df",
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -78L), spec = structure(list(
    cols = list(Run = structure(list(), class = c("collector_character",
    "collector")), `Assay Type` = structure(list(), class = c("collector_character",
    "collector")), AvgSpotLen = structure(list(), class = c("collector_double",
    "collector")), BioProject = structure(list(), class = c("collector_character",
    "collector")), BioSample = structure(list(), class = c("collector_character",
    "collector")), Cell_type = structure(list(), class = c("collector_character",
    "collector")), `Center Name` = structure(list(), class = c("collector_character",
    "collector")), Consent = structure(list(), class = c("collector_character",
    "collector")), `DATASTORE filetype` = structure(list(), class = c("collector_character",
    "collector")), `DATASTORE provider` = structure(list(), class = c("collector_character",
    "collector")), `DATASTORE region` = structure(list(), class = c("collector_character",
    "collector")), donorid = structure(list(), class = c("collector_character",
    "collector")), Experiment = structure(list(), class = c("collector_character",
    "collector")), GEO_Accession = structure(list(), class = c("collector_character",
    "collector")), Instrument = structure(list(), class = c("collector_character",
    "collector")), LibraryLayout = structure(list(), class = c("collector_character",
    "collector")), LibrarySelection = structure(list(), class = c("collector_character",
    "collector")), LibrarySource = structure(list(), class = c("collector_character",
    "collector")), MBases = structure(list(), class = c("collector_double",
    "collector")), MBytes = structure(list(), class = c("collector_double",
    "collector")), Organism = structure(list(), class = c("collector_character",
    "collector")), Platform = structure(list(), class = c("collector_character",
    "collector")), ReleaseDate = structure(list(format = ""), class = c("collector_datetime",
    "collector")), sample_acc = structure(list(), class = c("collector_character",
    "collector")), `Sample Name` = structure(list(), class = c("collector_character",
    "collector")), sample_type = structure(list(), class = c("collector_character",
    "collector")), source_name = structure(list(), class = c("collector_character",
    "collector")), `SRA Study` = structure(list(), class = c("collector_character",
    "collector"))), default = structure(list(), class = c("collector_guess",
    "collector")), skip = 1), class = "col_spec"))

```

I want to have the below output to be in my final dataframe which i want to use for renaming the original file.

```
 Run        source_name
   <chr>      <chr>
 1 SRR2920466 HSC1
 2 SRR2920477 HSC2
 3 SRR2920478 HSC3
 4 SRR2920506 HSC4
 5 SRR2920507 HSC5
 6 SRR2920532 HSC6
 7 SRR2920531 HSC7
 8 SRR2920469 CMP1
 9 SRR2920470 CMP2
10 SRR2920481 CMP3

```

Found a solution but I have to apply it to each string in the last column

```
i1 <- !sapply(df1, is.numeric)
df1[i1] <- lapply(df1[i1], function(x) factor(replace(as.character(x),
                                                      x == "HSC", paste0(x[x=="HSC"], seq_len(sum(x == "HSC"))))))

```

Right now I am able to match repeated string in the column and add increment.How to apply it for all the other string .

# Answer
We need to group by *source\_name* and add row id, here is using *data.table* package:

```
library(data.table)

# convert to data.table
setDT(df1)

# group by source_name, and add rowid
df1[, myNewColumn := paste0(source_name, rowid(source_name))]

df1
#            Run source_name myNewColumn
#  1: SRR2920466         HSC        HSC1
#  2: SRR2920477         HSC        HSC2
#  3: SRR2920478         HSC        HSC3
#  4: SRR2920506         HSC        HSC4
#  5: SRR2920507         HSC        HSC5
#  6: SRR2920532         HSC        HSC6
#  7: SRR2920531         HSC        HSC7
#  8: SRR2920469         CMP        CMP1
#  9: SRR2920470         CMP        CMP2
# 10: SRR2920481         CMP        CMP3
# ...

```

See here for more solutions.
> 3 votes

# Answer
EDIT: Wheres the OP accepted this answer, @zx8754's answer below, as well as the ones in the link they gave, are much more elegant and straightforward.

---

Here is a work-around with `table()`, `lapply()` and `seq_len()` (I have basically used your `lapply()` call on the output of a `table()` call):

```
my_counts <- table(df1$source_name)

> my_counts

Blast   CMP   GMP   HSC   LSC  Mono  pHSC
   23     8     7     7     8     6    19

df1 <- df1[order(df1$source_name),]

> head(df1)
          Run source_name
37 SRR2920546       Blast
38 SRR2920547       Blast
39 SRR2920548       Blast
40 SRR2920549       Blast
41 SRR2920550       Blast
42 SRR2920551       Blast

> names(unlist(lapply(my_counts, seq_len)))
 [1] "Blast1"  "Blast2"  "Blast3"  "Blast4"  "Blast5"  "Blast6"  "Blast7"  "Blast8"  "Blast9"  "Blast10" "Blast11" "Blast12" "Blast13"
[14] "Blast14" "Blast15" "Blast16" "Blast17" "Blast18" "Blast19" "Blast20" "Blast21" "Blast22" "Blast23" "CMP1"    "CMP2"    "CMP3"
[27] "CMP4"    "CMP5"    "CMP6"    "CMP7"    "CMP8"    "GMP1"    "GMP2"    "GMP3"    "GMP4"    "GMP5"    "GMP6"    "GMP7"    "HSC1"
[40] "HSC2"    "HSC3"    "HSC4"    "HSC5"    "HSC6"    "HSC7"    "LSC1"    "LSC2"    "LSC3"    "LSC4"    "LSC5"    "LSC6"    "LSC7"
[53] "LSC8"    "Mono1"   "Mono2"   "Mono3"   "Mono4"   "Mono5"   "Mono6"   "pHSC1"   "pHSC2"   "pHSC3"   "pHSC4"   "pHSC5"   "pHSC6"
[66] "pHSC7"   "pHSC8"   "pHSC9"   "pHSC10"  "pHSC11"  "pHSC12"  "pHSC13"  "pHSC14"  "pHSC15"  "pHSC16"  "pHSC17"  "pHSC18"  "pHSC19"

df1$new_col <- names(unlist(lapply(my_counts, seq_len)))

> head(df1)
          Run source_name new_col
37 SRR2920546       Blast  Blast1
38 SRR2920547       Blast  Blast2
39 SRR2920548       Blast  Blast3
40 SRR2920549       Blast  Blast4
41 SRR2920550       Blast  Blast5
42 SRR2920551       Blast  Blast6

```

You can store the initial row order somewhere to have resulting data frame in the same order the initial one.
> 1 votes

---
Tags: r
---


# Question
Title: How to generate conformations for a ligand for in silico screening?
I have multiple ligands (in pdb format) which I want to try and dock to various proteins.

But for each ligand I have only a single conformation (by which I meant the spatial relationship between it's atoms, not its position relative to the target protein).

How do I generate other possible conformations? Is there a generally used tool for that?

# Answer
> 1 votes
One option I found was Frog2

https://mobyle.rpbs.univ-paris-diderot.fr/cgi-bin/portal.py#forms::Frog2

https://github.com/tuffery/Frog2

From an input of a `mol2` or `sdf` file it creates a plethora of conformers, which can be outputted as a `pdb` file

# Answer
> 1 votes
## Open babel

The most commonly used tool is `open-babel` (obabel). It seems annoying that you have to specify the number of conformers, but this is because if you get something like retinol, the number of conformers is huge. The command from the documentation is:

```
obabel startingConformer.mol -O ga_conformers.sdf --conformer --nconf 30
   --score rmsd --writeconformers

```

An alternative argument is `--confab` which will pick each representative of a RMSD cutoff cluster.

Do note that 2D representations and un-protonated molecules will need extra command line arguments.

## RDKit

You can also generate conformers in RDkit, but I think this is beyond your requirements, but I thought it may be good to mention the essential compchem python library.

---
Tags: pdb
---


# Question
Title: Clustered column chart for single cell data
I want to create a column plot for my single cell analysis that has 2 sample types, a normal and PDAC. I want to show, the percentage of different cell types per sample type.

How do I go about determining the percent of each cell type in each condition and then go about plotting the column plot? Also, is it possible to determine the SD from the dataset and add it?

**dput(head(all@meta.data))**

```
structure(list(orig.ident = c("PDAC", "PDAC", "PDAC", "PDAC",
"PDAC", "PDAC"), nCount_RNA = c(7945, 7616, 7849, 5499, 853,
1039), nFeature_RNA = c(2497L, 2272L, 2303L, 2229L, 509L, 588L
), percent.mt = c(2.63058527375708, 3.7421218487395, 4.9433048796025,
0.490998363338789, 0.234466588511137, 0.192492781520693), RNA_snn_res.0.5 = structure(c(5L,
5L, 5L, 5L, 1L, 6L), .Label = c("0", "1", "2", "3", "4", "5",
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
"17", "18", "19"), class = "factor"), seurat_clusters = structure(c(5L,
5L, 5L, 5L, 1L, 6L), .Label = c("0", "1", "2", "3", "4", "5",
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
"17", "18", "19"), class = "factor"), S.Score = c(0.0171539298241524,
-0.00599651468836135, 0.011058574403656, -0.0424413740104151,
0.03549521163095, -0.0243347616753686), G2M.Score = c(-0.0591387992008088,
0.0253275642795205, -0.0512402869839816, -0.0239076248512967,
-0.0263164126515597, -0.0339086517371782), Phase = structure(c(3L,
2L, 3L, 1L, 3L, 1L), .Label = c("G1", "G2M", "S"), class = "factor"),
    old.ident = structure(c(1L, 1L, 1L, 1L, 1L, 3L), .Label = c("Fibroblast",
    "T Cell", "Endothelial", "Tumor", "Stellate", "Macrophage",
    "B Cell", "Mast Cell", "Acinar", "Endocrine", "Exocrine"), class = "factor")), row.names = c("PDAC_AAACCCAGTCGGCTAC-1",
"PDAC_AAACGAAGTCCAGGTC-1", "PDAC_AAAGAACAGCAAGTGC-1", "PDAC_AAAGAACAGGATGAGA-1",
"PDAC_AAAGTGACATCCGTTC-1", "PDAC_AACAACCAGGAAGTAG-1"), class = "data.frame")

```

# Answer
You can make use of `table()` and `prop.table()` to calculate the number and proportions of your cell types with respect to tissue of origin.

`prop.table()` requires a table as input so you can start with something like `prop.table(table(object@meta.data$cell_type, object@meta.data$tissue)`. This will create an objet of class `table`.

Next step is to create a data frame with `as.data.frame(your_table_object)` and then to convert it to "long format", for example with `tidyr::pivot_longer()`. You can feed these data in long format to `ggplot()`.
> 0 votes

# Answer
If I understand your question correctly and you have `dplyr` and `ggplot` installed, you could try something like (if a is your data.frame)

```
a %>%
dplyr::group_by(orig.ident, old.ident) %>%
dplyr::summarise(count  = n()) %>%
ggplot(aes(x=orig.ident, y=count)) + geom_bar(stat='identity', aes(group=old.ident), position='dodge')

```
> 0 votes

---
Tags: r, ggplot2
---


# Question
Title: Ensembl protein identifiers from different assemblies
In order to obtain orthologs for human genes, I am working with InParanoid version 7.0. This version of the database uses Ensembl Protein Identifiers (**ENSP**) from ENSEMBL version v54, which is based on **NCBI 36** assembly of the human genome.

For my project, I need to have the Ensembl Gene Identifiers (**ENSG**) of these proteins, but in the last version of ENSEMBL, which is based on **GRCh38** assembly. I need to use this version of InParanoid for consistency with previous analyzes.

Thus, my question is, **could I just use biomaRt to map ENSP (from NCBI36) to ENSG (actual version)** or this approach is inherently wrong?

Is it possible to "convert" ENSP identifiers from NCBI36 to GRCh38 and then do the map to ENSG?. I am aware of liftOver but I am not sure it would work here.

# Answer
Those IDs are elderly! Ensembl 54 was 2009!

I would recommend using BioMart combined with the ID history converter. The ID history converter will convert old IDs to new, and BioMart will convert ENSPs to ENSGs. You can either use the ID history converter with the ENSPs first, then the current BioMart to get the ENSGs. Or you can use Ensembl 54 BioMart to convert the old ENSPs to old ENSGs, then update the ENSGs with the ID history converter. I would probably do both and compare results.
> 4 votes

---
Tags: sequence-annotation, ensembl, liftover
---


# Question
Title: Bedtools wrongly indicates zero coverage
I have a BAM file produced using BWA-MEM and I need to calculate the coverage of each exon in this sequencing. I tried using this command:

`bedtools coverage -a exome_sorted.bed -b exome_alignment.bam -sorted`

And the result was the following:

```
chr1    12595   12802   ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242  .       ref|DDX11L1     0       0       207     0.0000000
chr1    13163   13658   ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242  .       ref|DDX11L1     0       0       495     0.0000000
chr1    14620   15015   ref|WASH7P,ref|NR_024540,ens|ENST00000538476,ens|ENST00000438504,ens|ENST00000488147,ens|ENST00000541675,ens|ENST00000423562       .       ref|WASH7P      0       0       395     0.0000000
...

```

That is, zero coverage for all exons. But there are reads mapping to the exons, as I can see by opening my BAM file on IGV and looking at some exonic regions:

I have already tried sorting the BED file using `sort -k1,1 -k2,2n` (as recommended on the Bedtools documentation) and using `bedtools sort` with the `-faidx` option, but neither option solved my issue.

if I skip the `-sorted` option, the memory usage becomes way too high.

I have also checked if the chromosome names are the same in the BAM file and in the BED file, using `samtools view -H`. Both use `chr1,chr2,...,chr22,chrX,chrY`, although BAM has `chrM` and stuff like `chrUn_gl000214` which are not on the BED.

Thanks!

EDIT:

@ATpoint pointed out that all visible reads in this region have MAPQ0. I looked at another exon that did not have this issue, but Bedtools still outputs zero coverage. IGV picture below:

# Answer
If this is your output :

```
chr1    12595   12802   ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242  .       ref|DDX11L1     0       0       207     0.0000000

```

Then I guess your bedfile looks like this (the first 6 fields) :

* chrom : chr1
* FeatureStart : 12595
* FeatureEnd : 12802
* name : ref|DDX11L1,ref|NR\_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242
* score : .
* strand : ref|DDX11L1

The issue I see here is that you are giving the wrong information to the strand field of your bedfile. It should be either `+`, either `-` or `.` if you don't want to consider the strand information.

Hope it helps !
> 2 votes

---
Tags: bedtools, coverage
---


# Question
Title: Is it possible for coronavirus or SARS to be synthetic?
I have heard several conspiracy theories regarding the origin of the new coronavirus, 2019-nCov. For example that the virus and/or SARS were produced in a laboratory or were some variant of Middle Eastern respiratory syndrome (MERS), shipped via laboratory workers.

I am well aware bioinformatics has debunked many conspiracy theories involving infectious diseases, an important one being polio vaccination programs in Africa were the origin of HIV, for example here. Likewise, numerous conspiracy theories on the man-made origins of HIV were similarly overturned using bioinformatic studies.

Do you have bioinformatics evidence that would debunk the current conspiracists about coronaviruses?

# Answer
> 24 votes
The scenarios are impossible and would be laughable if they were not so serious. The evidence is in the phylogenetic trees. Its a bit like a crime scene when the forensics team investigate. We've done enough crime-scenes often going to the site, collecting the pathogen, sequencing and then analysis - (usually neglected diseases) without any associated conspiracy theories.

The key technical issue is coronaviruses are **zoonoses** , pathogens spread to humans from animal reservoirs and phylogenetic tree really helps understand the how the virus is transmitted.

*Trees*

1. The key thing about all the trees are **bats**. Bats lineages are present at every single point of the betacoronavirus phylogeny (tree), both as paraphyletic and monophyletic lineages, one example is this tree of betacoronaviruses here. Meaning the nodes connecting the branches of the tree to the "master-branch", represent common ancestors and these were almost certainly bat-borne cornaviruses. This is especially true for SARS and - here bat-viruses are EVERYWHERE.
2. The tree here also shows that SARS arose on independently on two occassions, again surrounded by bat lineages and 2019-nCov has emerged separately at least once, again associated with bats.
3. Finally the tree below is a figure from BioRxiv Zhou et al (2020) "Discovery of a novel coronavirus associated with the recent pneumonia outbreak in humans and its potential bat origin" shows the 2019-nCov lineage is a direct descendent of a very closely virus isolated from a bat (RaTG13\*). This is a really conclusive finding BTW.

**Note**, I don't normally present inline images, but it is such a nice finding (hint to reviewers) and BioRxiv is open access.

**Conspiracy theory 1**: *laboratory made virus*

Literally it would require someone passaging a new virus, with unknown human pathogenicity, and independently introducing all the earlier passages enmass across bat populations of China. They would then hope each lineage becomes an indepedent virus population before, then introducing the virus to humans. Thus when field teams of scientists go around using mist nets to trap the bats, buy them from markets, isolate the virus and sequence it they would find a beautiful, array of natural variation in the bat populations leading up to the human epidemics, that perfectly matches vast numbers of other viral zoonoses. Moreover, this would have to have happen substancially prior to SARS and 2019-nCov, because the bat betacoronaviruses have been known about prior both epidemics, viz. its simply not feasible.

**Biological explanation** *General* Bats are a reservoir host to vast numbers of pathogens, particularly viruses, including many alphaviruses, flaviviruses, rabies virus and beleived to be important in ebolavirus (I don't know about this) and even important to several eukaryotic parasites. It makes sense, they are mammals, so evolutionary much closer to us than birds for example, with large dispersal potential and roost in 'overcrowded' areas enable rapid transmission between bats.

*Technical* The trees show bats are the common ancestor of betacoronaviruses in particular for the lineage leading into the emergence of 2019-nCov and SARS, this is seen in this tree, this one and the tree above. The obvious explanation is the virus circulates endemically in bats and has jumped into humans. For SARS the intermediate host, or possible "vector" was civet cats.

The theory and the observations fit into a seamless biological answer.

**Conspiracy theory 2**: *Middle Eastern connection*

I heard a very weird conspiracy theory attempting to connect MERS with 2019-nCov. The theory was elaborate and I don't think it is productive to describe here.

**Biological explanation** All the trees of betacoronaviruses show MERS was one of the earliest viruses to diverge and is very distant from 2019-nCov, to the extent the theory is completely implausible. The homology between these viruses is 50%, so its either MERS or 2019-nCov. Its more extreme than mixing up yellow fever virus (mortality 40-80%) with West Nile virus (mortality \<\<0.1%), the two viruses are completely different at every level.

**What about errors?** Phylogeneticists can spot it a mile off. There are tell-tale phylogenetic signatures we pick up, but also we do this to assess 'rare' genetic phenomina. There is nothing 'rare' about the coronaviruses. The only anomaly is variation in the poly-A tail and that is the natural variation from in vitro time-series experiments. Basically we've looked at enough virses/parasites through trees, that have no conspiracy theories at all (often neglected diseases), and understand how natural variation operates - so a phylogenecist can shift the wheat from the chaff without really thinking about it.

**Opinion** The conspiracy theories are deeply misplaced, and the only connection I can imagine is its China. However, the Chinese have loads of viruses, influenza in particular which causes major pandemics, but that is a consequence of their natural ecology (small-holder farming) allowing the virus to move between reservoir hosts. I've not visited small-hold farms in China, but I have in other parts of the world and when you see them, you get it. The pigs, chickens (ducks .. China), dogs, horses and humans all living within 10 meters of each other.

**Conclusion** Shipping large numbers of bats to market, bat soup, raw meat from arboreal (tree-living) mammals such as civets that are sympatric to bats. Then consider the classical epidemiology in light of the phylogenetic data, which is very consistent, a single picture emerges that coronavirus is one of many zoonoses which has managed to transmit between patients.

**Summary** The fundamental point is the bioinformatics fit into the classical epidemiology of a zoonose.

\*, **Note** The bat coronavirus RaTG13 predates the 2019-nCov outbreak by 7 years. It is not even clear whether the virus has been isolated, i.e. could just be a RNA sequence.

"They have found some 500 novel coronaviruses, about 50 of which fall relatively close to the SARS virus on the family tree, including RaTG13—it was fished out of a bat fecal sample they collected in 2013 from a cave in Moglang in Yunnan province."

Cohen, "Mining coronavirus genomes for clues to the outbreak’s origins" Feb, 2020 Science Magazine,

---
Tags: phylogenetics, phylogeny, ngs, covid-19
---


title: Seurat clusters

text:
# Question
Title: Seurat clusters
I'm trying to show cell cycle regression and input sample metadata to my umap.

I've been following Satija Lab's tutorials and have generated my umap but now want to specify which sample corresponds to each cluster, if one sample is driving a certain cluster, etc. I also, would like yo show cell cycle regression in each cluster. I following the cell cycle regression tutorial but that only tells me the board overview of the cells in each cycle. I would want to highlight it in the umap.

How would I go about each of these tasks?

# Answer
You can use the `group.by` argument of the `DimPlot()` function. `group.by` accepts column names present in the `meta.data` column of the Seurat objects.
> 1 votes

---
Tags: scrnaseq, seurat
---


# Question
Title: Why does this naive Python paired Fastq reader work?
I wrote some throwaway code for ingesting paired-end reads from a Fastq file into memory.

```
def ingest(stream1, stream2):
    while True:
        record1, record2 = list(), list()
        for _ in range(4):
            record1.append(next(stream1))
            record2.append(next(stream2))
        yield record1, record2

```

I had initially put the `for` loop in a `try/Except` block, anticipating I'd need to explicitly break out of the endless `while` loop. However, the code seems to work correctly even after I removed it. I even "corrupted" the final read in one of the paired files, and the `ingest` function performs reasonably (the last read pair is never `yield`ed).

My question is...**why does this code work**, since I never explicitly break out of the loop?

Here's an example of how the code is invoked.

```
>>> with open('seqR1.fastq', 'r') as fh1, open('seqR2.fastq', 'r') as fh2:
...     counter = 0
...     for r1, r2 in interleave(fh1, fh2):
...         counter += 1
>>> counter
5000
>>> r1
['@read5000/1\n', 'AGGTCTGGGGACAGGTGTTTTACAGGCTCCTTTAGGGGTGGCAAGTCTCCCGGAAACATGATTAGTGAAGACACGTTAACAGCGAAACTGACAGAGAACCTCAACATTCTCTGGTAGGAAAAGCCAGGATCCTGGCAGAGGAAGCAACGATTCGAGCAGGGAACTGTGCAGGACAGAGAACAATGCGGCCCCAATTTTGATCCAGGGGTAGTTAGCTTAGCTAAGGAGTTCCCCAGTCCCTAAGTGGAGAGGGGATCCAGAACCCTACAATTTATTATTATATCTTAGCGGATATATCTAT\n', '+\n', 'CCCCCGGGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG@GGGGGGGFAGGGCGGGGGFGGGGGGGGGCGGGGGDGGGGFGGGGGGGGGCGGGGGGGGGGGGGGGCGGGGGGGGGGGGGGGGG@G<GGGGGGGGGGGGFGGFGFF:DGGFFCGGGGCGGFEGGGGGGEGGGFB:GGGDG=FGGGG,GAGGGDGGFGGGCD3GFG+G?G*FGGG<GCGGGGC@CD<FFGF*:+G=FE<;GG,8C7GF=;9D=G<C5C66G6C>7)FD*6G7D)9;*C:)4D5+9>@G*=0):)*;1*504\n']
>>> r2
['@read5000/2\n', 'TGTGGGGTTCTGGATCCCCTCTCCACTTAGGGACTGGGGAACTCCTTAGCTAAGCTAATTACCCCTGGATCAAAATTGGGGCCGCATTGTTCTCTGTCCTGCACAGTTCCCTGCTCGAATCGTTGCTTCCTCTGCCAGGATCCTGGCTTTTCCTACCAGAGAATGTTGAGGTTCTCTGTCATTTTCGCTGCTAATGTGTCTTCACTAATCGTGTTTCCGGGAGACTTGCCACCCCTAAAGGAGCTTGTTGAACACTTGTCACCAGACCTAAGTCCTTAGGGATGAACCTGTGTCTCGGCAC\n', '+\n', 'CCCCBGC9FGCD-GGGFGGGEFEFF@CCGFGGGGFF+DG6,F69+FGG+FGF<FCC;GGEBFGGEDECFGGG6GDG7BFG9GCG9G9GDGGGG<GDFFFGEGC,<GC=GFG:,:A:EE+GD<:E,FGF=FCAGEFF++G,G,<EGGFGFFA>*FG9EGFCFCGC,B>43GFF@E@2>9ED@,568D=A;C+=5@E/F53D>,6F*9@+,E,C=*0+*@56<*9@:/F,1,:*)*14*0CC=*1)*0*,*)0;+)*))/**)**;0;//9)1+)10.()5/)070*)7*C**3/)(2*****\n']

```

# Answer
Your code only runs prior to Python 3.7, now it raises `RuntimeError: generator raised StopIteration`. This is essentially a bug fix described in PEP 479 -- Change StopIteration handling inside generators. You can enable the new behaviour from Python 3.5 via `from __future__ import generator_stop`.

Note, it's not related to file-handles per se, you can reproduce it with a smaller example with a simple iterator:

```
it = iter([1,2])

def f(it):
    while True:
        yield next(it)

print(list(f(it)))

```

You should not write code relying on the bug
> 4 votes

# Answer
EDIT: **Chris\_Rands's answer is more authoritative and clarifies that this behavior is essentially a bug.** (I had only tested my example on Python 3.6.8 and hadn't seen PEP 479.)

When `next` runs out of lines, it raises `StopIteration`, and apparently that stops the while loop similarly to calling `break` without you needing to catch the exception. I couldn't see a definitive source on that but you can see it in action if you catch and re-raise the exception:

```
def ingest2(stream1, stream2):
    while True:
        try:
            record1, record2 = list(), list()
            for _ in range(4):
                record1.append(next(stream1))
                record2.append(next(stream2))
            yield record1, record2
        except StopIteration:
            print("StopIteration raised")
            raise

```

You can also try to catch it outside the loop, but it doesn't get re-raised at that point. (If you try raising one yourself in a while loop, though, it *doesn't* catch it in the same way, so there must be some extra Python magic going on. I can't see any difference between the two exception objects so I'd be curious what's going on there.) Neat example. You could also make a pair of FASTQ parsers using Biopython's SeqIO.parse for more complex/permanent situations if needed.
> 2 votes

---
Tags: python, fastq
---


# Question
Title: Can open-source software be peer-reviewed and published?
My colleague and I have developed a software tool intend to release it open-source. This tool is specifically for tasks in bioinformatics but we think it would be helpful for the wider community. Our institution will permit us to release it provided we get appropriate credit.

Thus we wish to publish it in peer-review. Is peer-review publication of bioinformatics software available? If so what is required to publish it?

# Answer
> 13 votes
There are whole journals based primarily around publishing open source tools. The primary example of that is "Bioinformatics", where a lot of the open-source tools are published. We've also had luck publishing in the Nucleic Acids Research yearly special webservers issue, since we make Galaxy wrappers around our tools.

You can also publish on tools in venues like Nature Methods, but it's vastly more difficult to get in there if your paper is purely on a tool and it's not something game-changing like `salmon` or `kallisto`.

Regarding what's actually needed to publish a tool, you will typically need the following:

* Decent documentation. I ask for this as a reviewer and am asked about it when I submit things.
* Easy installation. No one wants to spend two hours trying to compile your tool, so make sure there's a conda package or at least a docker container for it.
* Test data. Without this your reviewers and users will likely not be happy since they won't be able to try it out on something small.

For the actual submission you'll tend to need comparisons to other tools. I'm generally not a big fan of this, but you will find that some reviewers will demand it (assuming there are other tools that do something similarish to yours).

# Answer
> 9 votes
Yes, there is even a journal dedicated to this: The Journal of Open Source Software

Here the software itself is peer reviewed and a decent documentation and tests are required. Since the submission and review process is based on GitHub, publishing there is straightforward if your code is on GitHub already.

Note: I am not affiliated with this journal, I just published two software packages there a while ago.

---
Tags: publishing
---


title: Random networks in PPi

text:
# Question
Title: Random networks in PPi
I am working with an in silico model to measure the distances between drugs targets and disease proteins of a disease. I would want to compare these distances with the distances of random networks but I don't know which is the best way to do the comparison. For instance, if I have 300 drug targets and 500 disease proteins mapped in the interactome and I want to measure the distance between them and then compare to random networks with the same features.

How do I generate these random networks? Choosing the same number of drug targets (300) and disease genes (500) randomly in my referral network? Choosing the same number of nodes (300 and 500) and the same number of links?

I'd appreciate your help!

# Answer
If your goal is only to account for degree distribution (i.e. Matteo Ferla's comment), a simulation approach you could use is:

1. measure the degree (number of links) of each protein in your PPI
2. find a large set of proteins that are equal or similar in degree to each of your proteins of interest (you may have to use a range of degrees to get a large enough set for comparison, ideally \>=1000). We'll call these "comparable proteins" for each of your proteins of interest.
3. For each protein of interest, make all of your measurements also for its set of comparable proteins
4. Compare your measurement (shortest distance?) for the protein of interest to that for its comparable proteins. You could compute Z-scores from the parameters of the comparable protein distribution, or just empirical Monte Carlo p-values from where it falls in the distribution.

This will give you at least a heuristic null distribution for each protein. It should be relatively easy to code up in whatever your favorite language is.

Hope that helps-
> 1 votes

---
Tags: networks
---


# Question
Title: Writing a python module to run Infernal on a list of fasta files
I am very new to using python3 and am hoping to use it to reduce the time I need to spend running Infernal to look for a variety of RNAs in a large collection of bacterial genomes. I am hoping to write a short module in python that will let me sequentially run infernal using a single covariance model on a list of genome files.

At the moment I am trying to use subprocess.run() or subprocess.Popen() to do so.

```

import subprocess
import os
import glob

genomes = list(glob.glob(os.path.join("/path/to/genomes/", "*.fa")))

for x in files:
    subprocess.Popen(args=["cmsearch", "/path/to/covariancemodel/", "genomes"], shell=False)

```

my issue is I do not know a way to have "genomes" passed as the components of the list within the command.

Is this a way off approach to trying to automate this? If so is there a better way to do it or is there a way to simply pass each file in the list to the cmsearch command?

# Answer
> 1 votes
I think what you want to do is:

```
import subprocess
import os
import glob

genomes = glob.glob(os.path.join("/path/to/genomes/", "*.fa"))

for gfile in genomes:
    subprocess.Popen(args=["cmsearch", "/path/to/covariancemodel/", gfile], shell=False)

```

each pass through the for loop will take a one of the genome files and pass it into the cmsearch command.

Note:

1) glob.glob returns a list so no need to cast it again - https://docs.python.org/2/library/glob.html

---
Tags: python, rna
---


title: renaming IDs in gff3 with BCBio.GFF

text:
# Question
Title: renaming IDs in gff3 with BCBio.GFF
I wrote a script which should changes IDs in a GFF3 file. Unfortunately, the below script has two problems.

1. It attaches the new ID to `Parent` which leads that `Parent` contains the old and new id. How is it possible to keep only the new one?
2. How is it possible to access the chromosome name?

This how to run the script `python renameIDgff3.py --gff3 braker_utr.gff3 --prefix ACTG --output braker_utr-newID.gff3`

```
#!/usr/bin/python3
import click
from BCBio.GFF import GFFExaminer
from BCBio import GFF

@click.command()
@click.option('--gff3', help="Provide GFF3 file", required=True)
@click.option('--prefix', help="e.g. ASSCTG", required=True)
@click.option('--output', help="Keep GFF3 file", required=True)
def run(gff3, prefix, output):
    print("Hello")
    with open(output, "w") as out_handle:
        for rec in GFF.parse(gff3):

            for count, feature in enumerate(rec.features):
                print("count", count)
                print(feature)
                print(feature.qualifiers.get("Name"))
                print(feature.sub_features)
                print(feature.sub_features[0].qualifiers.get("Name"))
                print("!!!!change")
                feature.qualifiers["ID"] = prefix + str(count).zfill(6)

                print(feature.sub_features[0].qualifiers["ID"])
                id_extension =  feature.sub_features[0].qualifiers["ID"][0].split('.')[1]
                feature.sub_features[0].qualifiers["ID"] = prefix + str(count).zfill(6) + '.' + id_extension
                print(feature.sub_features[0].qualifiers["Parent"])
                print("-----------")

            GFF.write([rec], out_handle)

if __name__ == '__main__':
    run()

```

Input file:

```
NbV1Ch08    AUGUSTUS    gene    7015    29794   0.01    -   .   ID=g1;
NbV1Ch08    AUGUSTUS    mRNA    7015    29794   0.01    -   .   ID=g1.t1;Parent=g1
NbV1Ch08    AUGUSTUS    transcription_end_site  7015    7015    .   -   .   Parent=g1.t1;
NbV1Ch08    AUGUSTUS    three_prime_utr 7015    8531    0.2 -   .   ID=g1.t1.3UTR1;Parent=g1.t1
NbV1Ch08    AUGUSTUS    exon    7015    8747    .   -   .   ID=g1.t1.exon1;Parent=g1.t1;
NbV1Ch08    AUGUSTUS    stop_codon  8532    8534    .   -   0   Parent=g1.t1;
NbV1Ch08    AUGUSTUS    CDS 8532    8747    0.31    -   0   ID=g1.t1.CDS1;Parent=g1.t1
NbV1Ch08    AUGUSTUS    intron  8748    9191    0.49    -   .   Parent=g1.t1;
NbV1Ch08    AUGUSTUS    CDS 9192    9342    0.66    -   1   ID=g1.t1.CDS2;Parent=g1.t1
NbV1Ch08    AUGUSTUS    exon    9192    9342    .   -   .   ID=g1.t1.exon2;Parent=g1.t1;
NbV1Ch08    AUGUSTUS    intron  9343    9915    0.58    -   .   Parent=g1.t1;
NbV1Ch08    AUGUSTUS    CDS 9916    10006   0.71    -   2   ID=g1.t1.CDS3;Parent=g1.t1
NbV1Ch08    AUGUSTUS    exon    9916    10006   .   -   .   ID=g1.t1.exon3;Parent=g1.t1;
NbV1Ch08    AUGUSTUS    intron  10007   10101   0.74    -   .   Parent=g1.t1;
NbV1Ch08    AUGUSTUS    CDS 10102   10201   0.78    -   0   ID=g1.t1.CDS4;Parent=g1.t1
NbV1Ch08    AUGUSTUS    exon    10102   10201   .   -   .   ID=g1.t1.exon4;Parent=g1.t1;
NbV1Ch08    AUGUSTUS    intron  10202   10712   0.8 -   .   Parent=g1.t1;
NbV1Ch08    AUGUSTUS    CDS 10713   11107   0.11    -   2   ID=g1.t1.CDS5;Parent=g1.t1
NbV1Ch08    AUGUSTUS    exon    10713   11107   .   -   .   ID=g1.t1.exon5;Parent=g1.t1;
NbV1Ch08    AUGUSTUS    intron  11108   11569   0.07    -   .   Parent=g1.t1;
NbV1Ch08    AUGUSTUS    CDS 11570   12151   0.09    -   2   ID=g1.t1.CDS6;Parent=g1.t1
NbV1Ch08    AUGUSTUS    exon    11570   12151   .   -   .   ID=g1.t1.exon6;Parent=g1.t1;
NbV1Ch08    AUGUSTUS    intron  12152   12588   0.34    -   .   Parent=g1.t1;
NbV1Ch08    AUGUSTUS    CDS 12589   12717   0.39    -   2   ID=g1.t1.CDS7;Parent=g1.t1
NbV1Ch08    AUGUSTUS    exon    12589   12717   .   -   .   ID=g1.t1.exon7;Parent=g1.t1;
NbV1Ch08    AUGUSTUS    intron  12718   12789   0.42    -   .   Parent=g1.t1;
NbV1Ch08    AUGUSTUS    CDS 12790   13075   0.39    -   0   ID=g1.t1.CDS8;Parent=g1.t1
NbV1Ch08    AUGUSTUS    exon    12790   13075   .   -   .   ID=g1.t1.exon8;Parent=g1.t1;
NbV1Ch08    AUGUSTUS    intron  13076   14832   0.51    -   .   Parent=g1.t1;
NbV1Ch08    AUGUSTUS    CDS 14833   15009   0.39    -   0   ID=g1.t1.CDS9;Parent=g1.t1
NbV1Ch08    AUGUSTUS    exon    14833   15009   .   -   .   ID=g1.t1.exon9;Parent=g1.t1;
NbV1Ch08    AUGUSTUS    intron  15010   15278   0.59    -   .   Parent=g1.t1;
NbV1Ch08    AUGUSTUS    CDS 15279   15415   0.56    -   2   ID=g1.t1.CDS10;Parent=g1.t1
NbV1Ch08    AUGUSTUS    exon    15279   15415   .   -   .   ID=g1.t1.exon10;Parent=g1.t1;
NbV1Ch08    AUGUSTUS    intron  15416   15487   0.58    -   .   Parent=g1.t1;
NbV1Ch08    AUGUSTUS    CDS 15488   15612   0.96    -   1   ID=g1.t1.CDS11;Parent=g1.t1
NbV1Ch08    AUGUSTUS    exon    15488   15612   .   -   .   ID=g1.t1.exon11;Parent=g1.t1;
NbV1Ch08    AUGUSTUS    intron  15613   15706   0.96    -   .   Parent=g1.t1;
NbV1Ch08    AUGUSTUS    CDS 15707   15957   0.98    -   0   ID=g1.t1.CDS12;Parent=g1.t1
NbV1Ch08    AUGUSTUS    exon    15707   15958   .   -   .   ID=g1.t1.exon12;Parent=g1.t1;
NbV1Ch08    AUGUSTUS    start_codon 15955   15957   .   -   0   Parent=g1.t1;
NbV1Ch08    AUGUSTUS    five_prime_utr  15958   15958   0.99    -   .   ID=g1.t1.5UTR1;Parent=g1.t1
NbV1Ch08    AUGUSTUS    five_prime_utr  27458   28250   0.37    -   .   ID=g1.t1.5UTR2;Parent=g1.t1
NbV1Ch08    AUGUSTUS    exon    27458   28250   .   -   .   ID=g1.t1.exon13;Parent=g1.t1;
NbV1Ch08    AUGUSTUS    five_prime_utr  29272   29794   0.08    -   .   ID=g1.t1.5UTR3;Parent=g1.t1
NbV1Ch08    AUGUSTUS    exon    29272   29794   .   -   .   ID=g1.t1.exon14;Parent=g1.t1;
NbV1Ch08    AUGUSTUS    transcription_start_site    29794   29794   .   -   .   Parent=g1.t1;

```

Output file:

```
##gff-version 3
##sequence-region NbV1Ch08 1 129222376
NbV1Ch08    AUGUSTUS    gene    7015    29794   0.01    -   .   ID=ACTG000000
NbV1Ch08    AUGUSTUS    mRNA    7015    29794   0.01    -   .   ID=ACTG000000.t1;Parent=g1,ACTG000000
NbV1Ch08    AUGUSTUS    transcription_end_site  7015    7015    .   -   .   Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    three_prime_utr 7015    8531    0.2 -   .   ID=g1.t1.3UTR1;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    exon    7015    8747    .   -   .   ID=g1.t1.exon1;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    stop_codon  8532    8534    .   -   0   Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    CDS 8532    8747    0.31    -   0   ID=g1.t1.CDS1;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    intron  8748    9191    0.49    -   .   Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    CDS 9192    9342    0.66    -   1   ID=g1.t1.CDS2;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    exon    9192    9342    .   -   .   ID=g1.t1.exon2;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    intron  9343    9915    0.58    -   .   Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    CDS 9916    10006   0.71    -   2   ID=g1.t1.CDS3;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    exon    9916    10006   .   -   .   ID=g1.t1.exon3;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    intron  10007   10101   0.74    -   .   Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    CDS 10102   10201   0.78    -   0   ID=g1.t1.CDS4;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    exon    10102   10201   .   -   .   ID=g1.t1.exon4;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    intron  10202   10712   0.8 -   .   Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    CDS 10713   11107   0.11    -   2   ID=g1.t1.CDS5;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    exon    10713   11107   .   -   .   ID=g1.t1.exon5;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    intron  11108   11569   0.07    -   .   Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    CDS 11570   12151   0.09    -   2   ID=g1.t1.CDS6;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    exon    11570   12151   .   -   .   ID=g1.t1.exon6;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    intron  12152   12588   0.34    -   .   Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    CDS 12589   12717   0.39    -   2   ID=g1.t1.CDS7;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    exon    12589   12717   .   -   .   ID=g1.t1.exon7;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    intron  12718   12789   0.42    -   .   Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    CDS 12790   13075   0.39    -   0   ID=g1.t1.CDS8;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    exon    12790   13075   .   -   .   ID=g1.t1.exon8;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    intron  13076   14832   0.51    -   .   Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    CDS 14833   15009   0.39    -   0   ID=g1.t1.CDS9;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    exon    14833   15009   .   -   .   ID=g1.t1.exon9;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    intron  15010   15278   0.59    -   .   Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    CDS 15279   15415   0.56    -   2   ID=g1.t1.CDS10;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    exon    15279   15415   .   -   .   ID=g1.t1.exon10;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    intron  15416   15487   0.58    -   .   Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    CDS 15488   15612   0.96    -   1   ID=g1.t1.CDS11;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    exon    15488   15612   .   -   .   ID=g1.t1.exon11;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    intron  15613   15706   0.96    -   .   Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    CDS 15707   15957   0.98    -   0   ID=g1.t1.CDS12;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    exon    15707   15958   .   -   .   ID=g1.t1.exon12;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    start_codon 15955   15957   .   -   0   Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    five_prime_utr  15958   15958   0.99    -   .   ID=g1.t1.5UTR1;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    five_prime_utr  27458   28250   0.37    -   .   ID=g1.t1.5UTR2;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    exon    27458   28250   .   -   .   ID=g1.t1.exon13;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    five_prime_utr  29272   29794   0.08    -   .   ID=g1.t1.5UTR3;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    exon    29272   29794   .   -   .   ID=g1.t1.exon14;Parent=g1.t1,ACTG000000.t1
NbV1Ch08    AUGUSTUS    transcription_start_site    29794   29794   .   -   .   Parent=g1.t1,ACTG000000.t1

```

Thank you in advance

# Answer
Removing the old id from 'Parent\`:

```
for subfeature in feature.sub_features:
    subfeature.qualifiers["Parent"].pop()

```

Chromosome name:

```
 for rec in GFF.parse(gff3):
            chrID = rec.id

```
> 1 votes

---
Tags: sequence-annotation, python, gff3
---


# Question
Title: How to specify input resfile for Rosetta?
I'm trying to run the following Rosetta command.

```
rosetta_scripts.linuxgccrelease -s data/binding/dataset/test/protein.pdb \
    -parser:protocol data/binding/test/scripts/pack.xml \
    -out:path:all data/binding/test/output/ \
    -resfile data/binding/test/output/resfile.txt \
    -out:prefix RosettaOutput_ \
    -nstruct 1 -ignore_zero_occupancy false -overwrite

```

However, I keep getting the error:

```
core.pack.task.operation.TaskOperations: {0} [ ERROR ] File:resfile.txt not found!

ERROR: Cannot open file resfile.txt
ERROR:: Exit from: /home/sean/src/rosetta_src_2019.35.60890_bundle/main/source/src/core/pack/task/operation/TaskOperations.cc line: 1645
Error: {0} [ ERROR ] ERROR: Exception caught by JobDistributor while trying to get pose from job 'RosettaOutput_2MWY_0001'
Error: {0} [ ERROR ] Treating failure as bad input; canceling similar jobs

[ ERROR ]: Caught exception:

File: /home/sean/src/rosetta_src_2019.35.60890_bundle/main/source/src/core/pack/task/operation/TaskOperations.cc:1645
[ ERROR ] UtilityExitException
ERROR: Cannot open file resfile.txt

AN INTERNAL ERROR HAS OCCURED. PLEASE SEE THE CONTENTS OF ROSETTA_CRASH.log FOR DETAILS.

protocols.jd2.FileSystemJobDistributor: {0} job failed, reporting bad input; other jobs of same input will be canceled: RosettaOutput_2MWY_0001
protocols.jd2.JobDistributor: {0} no more batches to process...
protocols.jd2.JobDistributor: {0} 1 jobs considered, 1 jobs attempted in 13 seconds

[ ERROR ]: Caught exception:

File: /home/sean/src/rosetta_src_2019.35.60890_bundle/main/source/src/protocols/jd2/JobDistributor.cc:327
1 jobs failed; check output for error messages

AN INTERNAL ERROR HAS OCCURED. PLEASE SEE THE CONTENTS OF ROSETTA_CRASH.log FOR DETAILS.

```

The error resolves itself if I move `resfile.txt` to the current working directory. Why can't `resfile.txt` be found?

# Answer
> 2 votes
The problem is my `pack.xml` specified a `resfile.txt` location which over-wrote my command-line option.

```
    <TASKOPERATIONS>
       <ReadResfile name="Design_Resfile" filename="resfile.txt"/>
    </TASKOPERATIONS>

```

If you remove the `filename` option in the XML file, you are then able to specify the path.

```
    <TASKOPERATIONS>
       <ReadResfile name="Design_Resfile"/>
    </TASKOPERATIONS>

```

---
Tags: protein-structure, computational-biochemistry, rosetta
---


# Question
Title: WGCNA: Problem with selecting soft threshold
I have 25 tumor samples with counts data. Initially, I filtered out low expressed genes and then converted counts to `varianceStabilizingTransformation` using `deseq2` package. With this data I started using `WGCNA` for co-expression network analysis.

For selecting the soft threshold I see very strange plot. R2 cutoff is 0.8 and I see that none of the scale free topology model fit is above that.

Here is the code I used:

`df` is a dataframe with genes as rows and 25 samples as columns with counts data.

```
library("DESeq2")
filtered.counts <- df[rowSums(df==0)<5, ]

U3 <- as.matrix(filtered.counts)
vsd <- vst(U3, blind=FALSE)

oed <- vsd

gene.names=rownames(oed)
trans.oed=t(oed)
dim(trans.oed)

n=16462;
datExpr=trans.oed[,1:n]
dim(datExpr)

SubGeneNames=gene.names[1:n]

library(WGCNA)
options(stringsAsFactors = FALSE);
allowWGCNAThreads()

powers = c(c(1:10), seq(from = 12, to=20, by=2));
sft=pickSoftThreshold(datExpr,dataIsExpr = TRUE,
                      powerVector = powers,corFnc = cor,
                      corOptions = list(use = 'p'),networkType = "unsigned")

   Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
1      1 0.000273  0.041          0.786 4140.00  4000.000 6730.0
2      2 0.428000 -1.130          0.852 1570.00  1400.000 3770.0
3      3 0.673000 -1.540          0.882  729.00   579.000 2410.0
4      4 0.737000 -1.720          0.891  383.00   268.000 1670.0
5      5 0.745000 -1.830          0.886  220.00   134.000 1210.0
6      6 0.704000 -1.990          0.860  134.00    71.700  909.0
7      7 0.737000 -1.980          0.890   85.90    40.300  701.0
8      8 0.742000 -2.020          0.903   57.30    23.500  551.0
9      9 0.733000 -2.090          0.917   39.40    14.200  441.0
10    10 0.750000 -2.080          0.934   27.80     8.810  357.0
11    12 0.770000 -2.080          0.952   14.80     3.670  243.0
12    14 0.775000 -2.090          0.951    8.43     1.660  171.0
13    16 0.397000 -2.820          0.459    5.06     0.801  124.0
14    18 0.409000 -2.770          0.474    3.18     0.406   91.4
15    20 0.421000 -2.720          0.490    2.06     0.215   68.8

# Plot the results
sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;

# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit, signed R^2",type="n", main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],labels=powers,cex=cex1,col="red");

```

Is this the right way? Which soft threshold power should I select?

Any help is appreciated. thanq.

# Answer
Don't worry about it too much. I would go with power 8 based on my general experience (also reflected in WGCNA FAQ) and on the mean connectivity around 50 and median around 20, which seems reasonable to me.
> 2 votes

---
Tags: r, rna-seq, networks, rna, wgcna
---


# Question
Title: Looking for a way to download all pdb files that detail a protein-drug complex
I am looking for a way to download all pdb files that detail a protein-drug complex.

While I could download the entirety of the pdb and search manually (or maybe write a script to discard all proteins that do not have an HETATM entry that details a non-water, non-salt molecule), that sounds effort intensive and error prone.

It seems reasonable that this subset already exist, but I was not able to find it.

Does such a list/database exists?

# Answer
Proving a negative is hard, but to my knowledge I do not know a database that has *all* the structure of *drug* bound protein. I can give some close matches.

First a few terms to be on the same page. A **ligand** is something that binds, which can be physiological, such as a **cofactor** or **substrate**, or not (especially man-made), such as a **drug**, i.e. a compound that has a physiological effect, or a **hit**, a small fragment from a fragment screen that will not bind with nanomolar affinity. Generally, ligands are small molecules, but may include oligopeptides and oligonucleotides. In the case of crystal structures, you also have **crystallisation additives**, which are not ligands, but are hetatoms in PDBs. Also, the definition of drug may vary. drugbank.ca has the following categories of drugs "approved", "experimental", "withdrawn" for pharmaceuticals, "illicit" for narcotics and "investigational" for probes and counts macromolecules as drugs.

* RCSB PDB: drug mapping table is a simple searchable table with ligands that are clinical drugs only —experimental compounds (=probes) will be absent.
* bindingdb.org contains many table, including a '3D' one. This table contains all the PDBs with ligands and their affinities. This contains crystallisation additives and natural ligands. Peter Curran's repo for the HotSpots algorithm (first written by Chris Radoux) has a blacklist, but you are bound to get stuff missed.
* bindingmoad.org is a database with all the protein-ligand structures from PDB up until 2018, but it includes natural ligands, like cofactors and substrates as above.
* https://klifs.vu-compmedchem.nl a well curated database of kinases, which has a lot of members and many ligands, but is not at all diverse. As a result it is a clean dataset often used as a testbed for novel methodologies.
* I mentioned in the comments the SGC Target Enabling Packages as a small subset of proteins with many hits each. But I realise that there is not a handy PDB list and the targets are not drugs, but probes.
> 2 votes

---
Tags: pdb, interactions
---


title: tp53 example in Momand's book

text:
# Question
Title: tp53 example in Momand's book
I am new to bioinformatics. While reading Momand's '*Concepts in Bioinformatics and Genomics*' , I understood clearly that each codon translates to a specific amino acid.

But then, he talks about the *TP53* gene, and then without explanation talks about P53... Is that a naming convention? It's not explicitly explained... Does that applies for any case? ie *TP90* gene encodes protein P90?

I don't understand the naming: correct me if wrong, but they named the gene after the protein's weight, which is 53 kilodaltons? (It confuses me with the location of the gene).

Lastly I do not understand: if there is a mutation in the DNA, the mutation shall appear in only one cell I guess, how then the others cells' DNA start to mutate? or is it that all DNA in all cells are mutated instantly?

# Answer
> 3 votes
This is possibly more of a molecular genetics question appropriate to biology.stackexchange.com, also note that it's usually better on stackexchange sites to restrict your post to a single question and do some preliminary research via google or your search engine of choice. But I don't see any harm in answering it here nonetheless.

The real answer is that naming conventions for genes/proteins are incredibly inconsistent and counterintuitive. You can't really draw any conclusions about a gene/protein based on its name alone. For example, some proteins such as p53 are named for biochemical features, others are named for phenotypes observed when their gene is mutated, such as for instance the gene *decapentaplegic* in fruit flies. Often they involve little jokes, such as the gene called *mothers against decapentaplegic* which represses the *decapentaplegic* mutant phenotype. So it's all a bit of a mess and people kind of assume that you'll follow along, which can be confusing.

It's a little rarer that gene names differ substantially from the names of the proteins they encode, but not unheard-of. This is indeed the case with your *TP53*/p53 example. See a short explainer here or a longer one here. What you'll notice is that there are in fact multiple different names for the same gene/protein. What you'll often see is that the shorter name (p53) is the one that people tend to use colloquially.

For your second question regarding how mutations are propagated, the short answer is that mutations are inherited from mother to daughter cells. Cells don't generally "infect" other cells with mutations (though in biology there are always exceptions). I think that what you are interested in is the difference between "somatic" and "germ-line" mutation.

Hope that helps.

---
Tags: gene
---


# Question
Title: Error in createMAE function: non-unique values when setting 'row.names' in TCGA LIHC Data
I am working with the `createMAE()` function of the ELMER Bioconductor package. While executing the `createMAE()` function using TCGA LIHC DNA Methylation 450K and RNASeq HTSeq - Counts data, I am getting this error:

```
Error in.rowNamesDF <-(x, value = value) :

duplicate 'row.names' are not allowed

In addition: Warning message:

non-unique values when setting 'row.names': 'TCGA-DD-AACA-02A', 'TCGA-DD-AACA-02B'

```

Can anyone tell me how to get rid of such error?

# Answer
> 1 votes
The error message is pretty clear: R requires unique row names on data frames and some of the "values" that you are trying to set as row names are non-unique.

You will need to provide a vector i) of length the same as your row number, ii) that is composed of non-unique values.

---
Tags: r, rna-seq, bioconductor, methylation
---