text:
Chaning string in rows
I have a too big data frame
> dim(anno_maf)
[1] 2635299     132
>

> head(anno_maf[1:3,1:4])
                     Tumor_Sample_Barcode Chromosome Start_Position End_Position
1: LP2000107-DNA_A01_vs_LP2000102-DNA_A01       chr5      137536800    137536803
2: LP2000330-DNA_A01_vs_LP2000316-DNA_A01      chr11       64012229     64012231
3: LP2000330-DNA_A01_vs_LP2000316-DNA_A01      chr17       41720917     41720917
>

> unique(anno_maf$Tumor_Sample_Barcode)
  [1] "LP2000107-DNA_A01_vs_LP2000102-DNA_A01" "LP2000330-DNA_A01_vs_LP2000316-DNA_A01" "LP2000332-DNA_A01_vs_LP2000323-DNA_A01"
  [4] "LP2000333-DNA_A01_vs_LP2000324-DNA_A01" "LP6005334-DNA_A02_vs_LP6005333-DNA_A02" "LP6005334-DNA_C01_vs_LP6005333-DNA_C01"

I want to remove everything thing from vs onward in anno_maf$Tumor_Sample_Barcode, for example leaving LP2000107-DNA_A01 from LP2000107-DNA_A01_vs_LP2000102-DNA_A01 but I don't know how to do that
Any help?

sub() is the way to go.
Your pattern is your regular expression: Check for "_vs_" followed by any character (provided with .) for an undefined length (provided with *).
my_text <- "LP2000107-DNA_A01_vs_LP2000102-DNA_A01"

sub(pattern = "_vs_.*",
    replacement = "",
    x = my_text)

[1] "LP2000107-DNA_A01"

Thank you, how I execute this on my data frame?

your_df$your_col <- sub(your_pattern, replacement_string, your_df$your_col)

For learners, it is always better to use parameters in the name = value format. I'd recommend sub(pattern = your_pattern, replacement = replacement_string, x = your_df$your_col), so everything is explicit and the order of the parameters doesn't matter.

@Ram RS, I agree and that is why the sub() bit is written that way in the original answer. And actually I try to follow Hadley Wickham's style guide, in which the name = value format is used from the second parameter on. sub(), however, unconventionally (at least for me) takes x as the third argument, hence I used the name = value pairs for all of the arguments used. For those who are interested, here is the link for the tidyverse style guide: https://style.tidyverse.org

You're right, sub, gsub, grep, etc have weird parameter ordering and hence highly benefit from this format.


text:
Convert FASTA to FASTQ with dummy quality scores
I have a FASTA file which I would like to convert into FASTQ format as the tool I want to use my data in requires it in FASTQ format. So dummy quality scores are fine.
Note: I am not using any functionality of the tool that will require the quality score.

One way of doing this is with two subcommands from the pyfastaq suite.
fastaq to_fake_qual in.fasta - | fastaq fasta_to_fastq in.fasta - out.fastq

The first tool, to_fake_qual, creates fake quality scores (default 40) for each base and the - sends that file (.qual) to stdout. The second tool, fasta_to_fastq, consumes both the original fasta and the quality scores coming from stdin and turns these into a fastq file.

A simple Biopython solution:
from Bio import SeqIO

for r in SeqIO.parse("myfile.fa", "fasta"):
    r.letter_annotations["solexa_quality"] = [40] * len(r)
    print(r.format("fastq"), end='')

Example:
$ cat << EOF > myfile.fa
>1
ATG
>2
ATGTAGA
EOF
$ python3 myscript.py
@1
ATG
+
III
@2
ATGTAGA
+
IIIIIII

I think the python package bioconvert does what you want if you don't provide qualities: https://github.com/bioconvert/bioconvert/blob/master/bioconvert/fasta2fastq.py
The following outputs lots of warnings for me, but does generate a fastq file test.fq from an existing fasta file test.fa (using dummy quality "I", and adding "None" to sequence headers that do not have a comment part in the original fasta):
bioconvert fasta2fastq test.fa test.fq

(This requires installing bioconvert using pip beforehand: pip install bioconvert.)

Here's a quick one-liner with no dependencies other than the UNIX core utils and Perl.
cat seqs.fasta | paste - - | perl -ne 'chomp; s/^>/@/; @v = split /\t/; printf("%s\n%s\n+\n%s\n", $v[0], $v[1], "B"x length($v[1]))' > seqs.fastq

Caveat: assumes your Fasta sequences are not wrapped, i.e. two lines per record.

I have just uploaded a Python tool to convert fasta to fastq: https://github.com/Mostafa-MR/A_tool_to_convert_fasta_to_fastq
Please note this is presented directly from a Jupyter notebook, i.e. it is a ipynb.

source: Stack Exchange

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

sort: votes


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

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

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

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

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

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

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

That gives me this figure

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

Just what I wanted, thank you!

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


text:
Stacked barplot for single cell analysis
I want to create a stacked bar-graph with different cell cycles for each cell type within each condition. I have uploaded the file for it. Don't know how to go about it as now I have another condition (cell cycle). I was able to do the bar graph comparing cell cycle per condition but now I want to compare cell type and cell cycle in each condition
output of dput(head(all.combined@metadata))
structure(list(orig.ident = c("treated", "treated", "treated",
    "treated", "treated", "treated"), nCount_RNA = c(1892, 307, 1348,
    3699, 4205, 4468), nFeature_RNA = c(960L, 243L, 765L, 1612L,
    1341L, 1644L), percent.mt = c(0.211416490486258, 1.62866449511401,
    4.45103857566766, 4.4065963773993, 0.0713436385255648, 3.87197851387645
    ), RNA_snn_res.0.5 = structure(c(11L, 11L, 5L, 6L, 11L, 13L), .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(11L, 11L, 5L, 6L, 11L, 13L), .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.476893835992198, -0.0200784617568548, -0.0335915198305002, -0.0247184276246385, 0.010785196602457, 0.0190008903712199), G2M.Score = c(0.204441469200986, 0.173804859670862, -0.0313235510969097, -0.0376796363661889, -0.0559526905696905, -0.0122031631356698), Phase = structure(c(3L, 2L, 1L, 1L, 3L, 3L), .Label = c("G1", "G2M", "S"), class = "factor"), old.ident = structure(c(7L,7L, 1L, 4L, 7L, 9L), .Label = c("Fibroblast", "T cell", "Macrophage", "Stellate", "Acinar", "Endothelial", "Tumor", "B cell", "Mast cell", "Ductal", "Islets of Langerhans"), class = "factor")), row.names = c("treated_AAACGCTAGCGGGTTA-1", "treated_AAAGGTAAGTACAGAT-1", "treated_AAAGTGAGTTTGATCG-1", "treated_AAATGGACAAAGTGTA-1",
    "treated_AACAAAGGTCGACTTA-1", "treated_AACAGGGTCCTAGCCT-1"), class = "data.frame")

output of dput(tail(all.combined@metadata))
structure(list(orig.ident = c("untreated", "untreated", "untreated",
"untreated", "untreated", "untreated"), nCount_RNA = c(901, 823,
1184, 1835, 1147, 1407), nFeature_RNA = c(482L, 479L, 649L, 1043L,
604L, 709L), percent.mt = c(1.77580466148724, 2.91616038882138,
4.22297297297297, 3.86920980926431, 2.0052310374891, 4.05117270788913
), RNA_snn_res.0.5 = structure(c(7L, 7L, 7L, 14L, 7L, 7L), .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(7L, 7L, 7L, 14L, 7L, 7L), .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.0320858200243315, 0.0304725660342869, 0.0215996091745327,
    0.0384166213301423, 0.144956251122548, -0.0242770509986111
    ), G2M.Score = c(0.0904224391544142, 0.050148242050667, -0.0178041670730754,
    -0.0112596867977946, -0.0519554524339088, -0.0136533184257381
    ), Phase = structure(c(2L, 2L, 3L, 3L, 3L, 1L), .Label = c("G1",
    "G2M", "S"), class = "factor"), old.ident = structure(c(5L,
    5L, 5L, 5L, 5L, 5L), .Label = c("Fibroblast", "T cell", "Macrophage",
    "Stellate", "Acinar", "Endothelial", "Tumor", "B cell", "Mast cell",
    "Ductal", "Islets of Langerhans"), class = "factor")), row.names = c("untreated_TTTGGTTGTCTAATCG-18",
"untreated_TTTGGTTTCCCGAGGT-18", "untreated_TTTGTTGAGAACTGAT-18",
"untreated_TTTGTTGAGCTCGGCT-18", "untreated_TTTGTTGAGTGCCTCG-18",
"untreated_TTTGTTGCACGGTGCT-18"), class = "data.frame")

Been using this code to generate the previous graph.
ggplot(CC, aes(x = Condition, y = Percent, fill = Cell_Cycle))+
  geom_bar(stat = "identity")+
  geom_text(aes(label = paste(round(Percent,2),"%")), position = position_stack(vjust =  0.5))

What is DF ? Can you provide an image of the graph when you are comparing cell cycle in each condition ? Your example "all.combined@metadata" contains only "untreated" conditions, maybe you should update your example to provide few lines of treated conditions too.

geom_bar(stat = "identity") is rarely correct, as bar plots should be used for count-level data only. If you want to use a value, then use geom_col() instead - this is explained in the help for geom_bar.

You don't have to post all the metadata file, but few lines of each condition "treated / untreated" could be helpful.

@dc37 I posted the head and tail of the dataset.

Based on your example, you can use count for each phase of cell cycle for each cell types in each condition. Here, I'm using dplyr package to do that but most likely, you could have the same output using various methods:
DF <- rbind(dfhead, dftail)

library(dplyr)
DF_Count <- DF %>%group_by(orig.ident,Phase,old.ident) %>%
  count() %>%
  ungroup() %>%
  group_by(orig.ident,old.ident) %>%
  mutate(Freq = n/sum(n)*100)

# A tibble: 8 x 5
# Groups:   orig.ident, old.ident [5]
  orig.ident Phase old.ident      n  Freq
  <chr>      <fct> <fct>      <int> <dbl>
1 treated    G1    Fibroblast     1 100
2 treated    G1    Stellate       1 100
3 treated    G2M   Tumor          1  33.3
4 treated    S     Tumor          2  66.7
5 treated    S     Mast cell      1 100
6 untreated  G1    Acinar         1  16.7
7 untreated  G2M   Acinar         2  33.3
8 untreated  S     Acinar         3  50

As you can see, DF_Count has the frequency for each cell types of each phase of the cell cycle in function of the condition. We can use DF_Count to get the following plot. Using facet_wrap, you can create two panels based on the condition column and thus represent the cell cycle of each cell types in function of the treatment condition:
ggplot(DF_Count, aes(x = old.ident, y = Freq, fill = Phase))+
  geom_col()+
  geom_text(aes(label = paste(round(Freq, 2),"%")),position = position_stack(vjust = 0.5))+
  facet_wrap(~orig.ident)

Does it look what you are expecting ?
NB: DF here is the addition of your head and tail of your all.combined dataframe.

How about the other cell types? Can I just make DF <- all.combined@metadata? That way all the cell types are included?

It should work. Here the graph represent only the portion of data you have provided. If you think this answer was useful, do not hesitate to vote or validate it. Learn more here: https://bioinformatics.stackexchange.com/help/someone-answers


text:
Is it possible to obtain RefSeq transcripts from Ensembl database for GRCh37 to match Ensembl Genes?
I know that with GRCh38, it is possible to extract the matching RefSeq transcripts for Ensembl genes/transcripts/proteins using Biomart from Ensembl.
However, for GRCh37, this option does not exists. Using the following MySQL command returns Ensembl transcript IDs for matching RefSeq transcripts:
mysql -u anonymous -h ensembldb.ensembl.org -A -P 3337 -D homo_sapiens_core_99_37 -e "SELECT transcript.stable_id, xref.display_label FROM transcript, object_xref, xref,external_db WHERE transcript.transcript_id = object_xref.ensembl_id AND object_xref.ensembl_object_type = 'Transcript' AND object_xref.xref_id = xref.xref_id AND xref.external_db_id = external_db.external_db_id AND external_db.db_name = 'RefSeq_mRNA';"

Also, I was trying to play around with the databases schema but I am not able to find any solution to problem.
Ultimately, I wish to obtain a table like the following:
Gene stable ID      Transcript stable ID    Protein stable ID   RefSeq match transcript
ENSG00000120329     ENST00000239451         ENSP00000239451     NM_031947.4
ENSG00000164219     ENST00000419445         ENSP00000404676     NM_005023.4
ENSG00000150594     ENST00000280155         ENSP00000280155     NM_000681.4
ENSG00000113578     ENST00000337706         ENSP00000338548     NM_000800.5

Thank you in advance.

The RefSeq match option in BioMart is from the Matched Annotation from NCBI and EBI (MANE) collaboration between RefSeq and Ensembl. It has only been calculated for the up-to-date gene annotation on GRCh38 so cannot be obtained on GRCh37. You can get mapping from Ensembl to RefSeq transcripts through BioMart as RefSeq mRNA ID (refseq_mrna in R) but this is not a perfect match like the MANE, it is a mapping based on sequence similarity and similar genomic location, and there can be mismatches between them.

Thank you @Emily_Ensembl for your answer. I will continue discussing this topic here. Is it possible to obtain also the RefSeq mRNA ID with the versions using biomaRt ?

No, we don't have the versions on the old verion of the database.

Try the GRCh37 version of BioMart, and the External reference "RefSeq mRNA ID" / "RefSeq ncRNA ID":
https://grch37.ensembl.org/biomart
From the Perl code provided by BioMart, it looks like the field name should be "refseq_mrna" or "refseq_ncrna". I'm not sure how that translates to the SQL code.
I feel compelled to add a caution to this, though. GRCh37 is a very old build, having been replaced by GRCh38 in December 2013. If possible, please find a way to use GRCh38 for your analysis instead, and don't create any new software that depends on GRCh37. For example, gnomAD v3 supports GRCh38 (see bioconductor, or here).

Thank you for your comment. Unfortunately, I have to use GRCh37 since gnomAD were annotated with this old version using VEP from Ensembl. So, doing a lift-over or annotating again would lose any supports that the paper had at that time.


text:
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.

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?
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

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


text:
Making a boxplot or violinplot for several dataframe with different number of samples
I have 3 data frames for treatment of 3 different drugs and in each of them I have the number of mutation types like insertion, deletion, SNP and total of mutations for each patients. In each group I have different number of patients like below
> head(dat1)
                patient  DEL  INS   SNP MNP total
1:    LP6008337-DNA_H06  927  773 40756   0 42456
2:    LP6008334-DNA_D02 1049  799 31009   0 32857
> dim(dat1)
[1] 21  6
>

 head(dat2)
                Patient   DEL  INS   SNP MNP total
1:    LP6008031-DNA_E01 13552 3374 62105   0 79031
2:    LP6005500-DNA_G01   539  500 43451   0 44490
> dim(dat2)
[1] 33  6
>

> head(dat3)
                Patient   DEL   INS    SNP MNP  total
1:    LP6005935-DNA_F03 39168 16739  58095   0 114002
2:    LP6008269-DNA_D08   849   910 103501  11 105271

> dim(dat3)
[1] 106   6

I want to show if the number of each mutational category and total is different for each drug by boxplot, violinplot or any type of visualization
For example this picture is very nice
Can you help me?

To help us help you, please make your data available with dput() next time.
To achieve your goal with ggplot2, you would need all of your data in one data frame and in the "long format".
dat1 <- "
Patient DEL INS SNP MNP total
LP6008337-DNA_H06 927 773 40756 0 42456
LP6008334-DNA_D02 1049 799 31009 0 32857
"

dat2 <- "
Patient DEL INS SNP MNP total
LP6008031-DNA_E01 13552 3374 62105 0 79031
LP6005500-DNA_G01 539 500 43451 0 44490
"

dat3 <- "
Patient DEL INS SNP MNP total
LP6005935-DNA_F03 39168 16739 58095 0 114002
LP6008269-DNA_D08 849 910 103501 11 105271
"

# read data and add a column indicating data frame of origin (drug in this case)
dat1 <- readr::read_delim(dat1, delim = " ")
dat1$drug <- "drug1"
dat2 <- readr::read_delim(dat2, delim = " ")
dat2$drug <- "drug2"
dat3 <- readr::read_delim(dat3, delim = " ")
dat3$drug <- "drug3"

# bind rows
dat <- rbind(dat1, dat2, dat3)

> dat
# A tibble: 6 x 7
  Patient             DEL   INS    SNP   MNP  total drug
  <chr>             <dbl> <dbl>  <dbl> <dbl>  <dbl> <chr>
1 LP6008337-DNA_H06   927   773  40756     0  42456 drug1
2 LP6008334-DNA_D02  1049   799  31009     0  32857 drug1
3 LP6008031-DNA_E01 13552  3374  62105     0  79031 drug2
4 LP6005500-DNA_G01   539   500  43451     0  44490 drug2
5 LP6005935-DNA_F03 39168 16739  58095     0 114002 drug3
6 LP6008269-DNA_D08   849   910 103501    11 105271 drug3

# convert data to long format, ggplot() expects that
dat_long <- tidyr::pivot_longer(dat, cols = 2:6)

> dat_long
# A tibble: 30 x 4
   Patient           drug  name  value
   <chr>             <chr> <chr> <dbl>
 1 LP6008337-DNA_H06 drug1 DEL     927
 2 LP6008337-DNA_H06 drug1 INS     773
 3 LP6008337-DNA_H06 drug1 SNP   40756
 4 LP6008337-DNA_H06 drug1 MNP       0
 5 LP6008337-DNA_H06 drug1 total 42456
 6 LP6008334-DNA_D02 drug1 DEL    1049

ggplot(dat_long,
                aes(x = name, y = value)) +
  geom_boxplot() +
  facet_wrap(~drug)

EDIT: setting scales argument of facet_wrap() to "free" could be helpful
ggplot(dat_long,
       aes(x = name, y = value)) +
  geom_boxplot() +
  facet_wrap(~drug, scales = "free")

Sorry @haci how we show if the groups are significantly different in terms of each feature let's say the number of SNP? I have made one by geom_signif but it compares features internally rather than between groups

@Angel, this should be a separate question but in general people use ANOVA to test the difference of three or more groups. You will have to make sure that your data meet its assumptions though.

Thank you by var.test function I compared them although I am not sure if t test is suitable

If you use a violin plot instead and ggplot2 to make it then:
geom_violin(...not relevant stuff..., scale="count")

will scale each violin plot according to the number of actual entries in it. I often find this nice so small groups look super skinny and people don't focus on random changes in them so much.

Sorry what you mean by small groups look super skinny in my 3 data types?

I didn't say anything about your groups, reread what I wrote.


text:
How do I get the compartment Id and volume from SBML using libSBML Python?
I searched the Python documentation for libSBML but couldn't find a staightforward answer to this question. For a given SBML model, how can I use Python libSBML to retrieve the values of the id and volume attributes for a given compartment?

There aren't any examples in the python docs on how to do this. The best bet is to search the C/C++ docs which are more extensive. There is a python example in the docs which shows how to retrieve the number of compartments in a model but it doesn't go further. The code below shows how you can retrieve the information requested, assuming you have a model called testmodel.xml:
import libsbml

document = libsbml.readSBML("testmodel.xml")

if document.getNumErrors() > 0:
   print("Encountered the following SBML errors:" + "\n")
   document.printErrors()
else:
   model = document.getModel()

   nCompartments = model.getNumCompartments()

   print ("Number of compartments: ", nCompartments)
   for i in range (nCompartments):
       comp = model.getCompartment(i)

       print ("Id = ", comp.getId())
       print ("Volume = ", comp.getVolume())


text:
Calculating mean accross rows with repeated entries in R
I have this expression matrix
> head(dat[1:8,1:4])
  Gene.Symbol        p1        p2         p3
1       STAT1     7.6159     7.6784     7.4430
2       STAT1     7.4435     7.2850     7.1250
3       STAT1     7.8925     7.9269     7.4956
4       STAT1     7.5867     7.9285     7.6741
5       GAPDH     8.4861     8.3189     9.0913
6       GAPDH     9.0206     8.0813     9.2360
> dim(dat)
[1] 58324   276
>

As you are seeing, I have repeated genes, so how I take average over expression values for redundant genes finishing with one unique expression value for each gene for each sample in columns?

Using group_by from dplyr
You can use group_by function from dplyr to calculate the mean of each gene symbol in each samples:
library(dplyr)
dat %>%
  group_by(Gene.Symbol) %>%
  summarise_at(vars(contains("p")), funs(mean(., na.rm = TRUE)))

# A tibble: 2 x 4
  Gene.Symbol    p1    p2    p3
  <chr>       <dbl> <dbl> <dbl>
1 GAPDH        8.75  8.20  9.16
2 STAT1        7.63  7.70  7.43

Using data.table
Using data.table to speed up the calculation, you can do:
library(data.table)
setDT(dat)
dat[, lapply(.SD, mean), by = Gene.Symbol]

   Gene.Symbol      p1     p2       p3
1:       STAT1 7.63465 7.7047 7.434425
2:       GAPDH 8.75335 8.2001 9.163650
```


text:
stacked bargraph with error bars
I'm attempting to plot a stacked barplot with ggplot2 with this code:
ggplot(CC, aes(x = Condition, y = Percent, fill = Cell_Cycle))+
  geom_bar(stat = "identity")+
  geom_text(aes(label = paste(round(Percent,2),"%")), position = position_stack(vjust =  0.5))

Question: How do I go about adding error bars to the graph?
I tried this code:
ggplot(DF, aes(x = Condition, y = Percent, fill = Cell_Type))+
     geom_bar(stat = "identity")+
     geom_text(aes(label = paste(round(Percent,2),"%")), position = position_stack(vjust =  0.5) + geom_errorbar(aes(ymin=Condition, ymax=Condition), width = .2))

But got this error:

Error: Cannot add ggproto objects together.

output of dput(DF)
structure(list(Cell_Type = structure(c(1L, 2L, 3L, 4L, 5L, 6L,
7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L), .Label = c("Fibroblast", "T cell",
"Macrophage", "Tumor", "Islets of Langerhans", "Endothelial",
"B cell"), class = "factor"), Condition = structure(c(1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("treated",
"untreated"), class = "factor"), Freq = c(6051L, 1892L, 1133L,
657L, 116L, 868L, 832L, 5331L, 3757L, 1802L, 835L, 287L, 704L,
256L), Percent = c(52.3941466793662, 16.3823707680319, 9.8103731924842,
5.68880422547407, 1.00441596675037, 7.51580223395965, 7.20408693393367,
41.0962072155412, 28.9623805118717, 13.8914585260561, 6.43694110391613,
2.21245760098674, 5.42707369719396, 1.97348134443417)), class = c("grouped_df",
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -14L), groups = structure(list(
    Condition = structure(1:2, .Label = c("treated", "untreated"
    ), class = "factor"), .rows = list(1:7, 8:14)), row.names = c(NA,
-2L), class = c("tbl_df", "tbl", "data.frame"), .drop = TRUE))

What are len and sd ? can you add a reproducible example of DF, len,sd  ? By curiosity, does your question is related to this one https://bioinformatics.stackexchange.com/questions/11197/stacked-barplot-for-single-cell-analysis?

I didn't change len and sd when I pasted it. I have updated it to Condition and added the output of DF. And yes, it is related to that question. Wasn't sure whether to ask it on the same question.

In your exemple, there is no mention of len and sd, so it is not clear what are these values and how you want to calculate them. Also, does the related question is yours ? because you don't have the same username ?

len and sd are not part of my example. When I pasted the code to the question, I didn't change it. Yes, it is my question, as for the username, I just created a new account as I couldn't log into the other.

Ok, so your question is about to calculate len and sd ? because you don't have these values, right ? For your account, you can see this link: https://bioinformatics.stackexchange.com/help/merging-accounts. Also, when you are registered on one site of StackExchange, you can directly create a new account on another StackExchange website by clicking on topright join Community button, it will automatically link your account on all sites you are registered of.

You are getting Error: Cannot add ggproto objects together. because you lack a ) and geom_text() is not closed. And once this is fixed, you would get Error: Discrete value supplied to continuous scale because geom_errorbar() is for intervals but your Condition variable is a factor.

@dc37 No, my question is how to do a stacked bar graph with the error bars.

@haci How do I change it to an interval?

@user6810 Well, you have to have some min and max values to plot. Depends on why you need these error bars.

I want to indicate the error or uncertainty of each cell type.

If you want to indicate the error or the uncertainty for each cell type, you need to replicate your experiments, ccount each cell type in each replicate, calculate the average percentage of each cell type in each condition and the associated standard deviation.
To answer your question about how to plot a stacked bargraph with error bar, first this option is not integrated into ggplot2 (see this discussion from the creator of ggplot2: https://github.com/tidyverse/ggplot2/issues/1079)
So, you need to calculate the position of each error bar by hand. You can do that by using dplyr. Here, I set the same error bar for everyone as you don't know what is the standard deviation in your data but at least you can observe the structure of the dataframe:
library(dplyr)
library(ggplot2)
DF %>% mutate(SD = 2) %>%
  mutate(Cell_Type = factor(Cell_Type, levels = c("B cell","Endothelial","Islets of Langerhans","Tumor","Macrophage","T cell","Fibroblast"))) %>%
  group_by(Condition) %>%
  mutate(SDPos = cumsum(Percent)) %>%
  ggplot(aes(x = Condition, y = Percent, fill = Cell_Type))+
  geom_bar(stat = "identity")+
  geom_text(aes(label = paste(round(Percent,2),"%")), position = position_stack(vjust =  0.5))+
  geom_errorbar(aes(ymin = SDPos-SD, ymax = SDPos+SD), width = 0.3, position = "identity")

As you can see, it's pretty ugly... so, if you need to add error bar, I will rather advised you to do a dodged barchart as the following:
library(dplyr)
library(ggplot2)
DF %>% mutate(SD = 2) %>%
  ggplot(aes(x = Cell_Type, y = Percent, fill = Condition))+
  geom_bar(stat = "identity", position = position_dodge())+
  geom_text(aes(label = paste(round(Percent,2),"%")), position = position_dodge(0.9), vjust =-2, size = 3)+
  geom_errorbar(aes(ymin = Percent-SD, ymax = Percent+SD), width = 0.3, position = position_dodge(0.9))+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  scale_y_continuous(limits = c(-2,60))

To my opinion, it makes easier the comparison of the distribution of each cell type between conditions.
Hope that it answers your question.

For this scenario, I would advise against a barchart and would use a violin plot or a boxplot instead. The latter two would enable to assess the number of replicates within a group as well as the spread of the replicates.

I agree if you are able to have more than 3 replicates. Otherwise, the spread observed with boxplot or violin plot won't mean a lot. But as it seems that individual values are the count of different cell types from single cell RNA seq, I doubt you can have a lot of replicates. That's why I advised for barchart. But it is also possible to plot individual values and the mean using geom_pointrange.

These numbers and graphs will not mean anything anyway if one does not have many replicates, especially if the data are from human samples and not from a model organism with a rather "clean" genetic background. Moreover, it is even more important to use box/violin plots when the number of replicates are small so that the reader is not mislead about the lack of power in the analyses. I guess, at least hope, that we will not see these plots in the future anyway: https://twitter.com/rafalab/status/1098973538260840449?lang=en

Totally agree ;) and personally, I won't publish those kind of results (at least in this format). But it's well beyond the scope of the question :D

These samples are derived from human samples, and I have replicates for each sample. A total of 32 samples, 16 from each condition. What do you all propose I should use to plot? I'm new at bioinformatics, and assumed a bargraph would be ideal. From reading the comments and answer, I don't have SD for each. I will add that.

If you have 16 for each conditions, you should associate them in a single dataframe and plot a boxplot or a violin plot as advised by @haci.


text:
Interpreting MultiQC plot of GATK BaseRecalibration data?
I'm using GATK4 to build a somatic variant calling pipeline. The pipeline uses MultiQC to aggregate quality control data, and one of the QC measures reported is base quality score recalibration from GATK4's BaseRecalibraton. The description in MultiQC says

This plot shows the distribution of base quality scores in each sample before and after base quality score recalibration (BQSR). Applying BQSR should broaden the distribution of base quality scores.

The black line is for my tumor sample and blue line is my normal sample. I don't understand how the before and after distributions are both shown in the plot as I only see 1 distribution per sample. Am I missing something, or should there be a "Post-calibration Count" plot that isn't shown? If nothing is missing, how should I be interpreting the plot?

Did you actually run the ApplyBQSR tool?

I ran both BaseRecalibrator and ApplyBQSR and the BAM headers show ApplyBQSR has been run. Is the post calibration data only stored in the BAM files or is it also in the tables output by BaseRecalibrator? If the post calibration data is only in the BAM files, I'm not sure how MultiQC is supposed to know the data exists...

There may be an option to have GATK output a metrics file.

Turns out two recalibration tables are needed to show pre-recalibration and post-recalibration plots. I previously only had generated the pre-calibration table by running BaseRecalibrator on the BAM that hadn't undergone ApplyBQSR to adjust base scores. Here are the steps to show pre and post recalibration plots in MultiQC:

Run BaseRecalibrator on your initial input BAM. This will create a recalibration table with the original base scores.
Run ApplyBQSR with the recalibration table from step 1 to create a new BAM with the adjusted base quality scores.
Run BaseRecalibrator on the adjusted BAM generated in step 2 to create a new recalibration table with the adjusted base quality scores.
Use both recalibration tables as inputs to MultiQC. This should create both plots.

hey did you method work ? I have done your steps and it does not create both plots. I believe we need to create a csv file report from gatk AnalyzeCovariates. Will try that now

@alex what part does not work? Does MultiQC not create both plots even though you give it both tables as input?

yes correct. the tables were in different directories and didn't work. I then tried to put the tables under the same directory, did not work either. then i noticed both tables had the same 'sample name'. I changed the names of the tables so multiqc would detect 2 'samples' and correctly use the tables as before and after. It also did not work. strangely their old test data for BQSR on github works. things have probaly changed since then and they did not update the code

Not sure what the issue is then. It has been awhile since I used that MultiQC function and eventually dropped those plots from my pipeline. I would ask the creator of MultiQC, Phil Ewels for help on his Gitter. There's instructions for how to contact him Github.


text:
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!

Why not gain these parameters from the distribution of real PPI networks? Say, https://www.ebi.ac.uk/intact/.

I have calculated these parameters from the real PPi but I want to compare with the random expectation to see if my datasets are significantly close between them...

A worry here I reckon is the distribution. I mean the model that you choose, say equiprobable vs. some form of power law.  In the most simplest system, a mix of ideal gases the collision frequencies will depend on the concentrations, so will be equiprobable only if the concentrations are equal. In reality, protein have different concentrations and interactions strengths. The distribution of these will strongly depend on the dataset —I don't suggest a good model that accounts for size etc. but just a fit to a gamma distribution.

Matteo Ferla, thanks for the comment but we don't have information about the strengths of the interactions and protein concentration. We only have an undirected and unweighted network with two different groups of proteins and I'm new in this topic so I don't understand what you said about fitting to a gamma distribution. :( I have seen that in most cases they compare their networks with random networks with the same number of nodes and/or links randomly choosen.

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

measure the degree (number of links) of each protein in your PPI
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.
For each protein of interest, make all of your measurements also for its set of comparable proteins
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-


text:
Why does the SARS-Cov2 coronavirus genome end in aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa (33 a's)?
The SARS-Cov2 coronavirus's genome was released, and is now available on Genbank.  Looking at it...

    1 attaaaggtt tataccttcc caggtaacaa accaaccaac tttcgatctc ttgtagatct
   61 gttctctaaa cgaactttaa aatctgtgtg gctgtcactc ggctgcatgc ttagtgcact
  121 cacgcagtat aattaataac taattactgt cgttgacagg acacgagtaa ctcgtctatc
  ...
29761 acagtgaaca atgctaggga gagctgccta tatggaagag ccctaatgtg taaaattaat
29821 tttagtagtg ctatccccat gtgattttaa tagcttctta ggagaatgac aaaaaaaaaa
29881 aaaaaaaaaa aaaaaaaaaa aaa

Wuhan seafood market pneumonia virus isolate Wuhan-Hu-1, complete genome, Genbank

Geeze, that's a lot of a nucleotides---I don't think that's just random.  I would guess that it's either an artifact of the sequencing process, or there is some underlying biological reason.
Question: Why does the SARS-Cov2 coronavirus genome end in 33 a's?

Not a biologist, but it seems awfully similar to the mechanism in of a NOP sled, just a naturally occurring one from nature. https://en.wikipedia.org/wiki/NOP_slide

29821 - 29881 (aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa)

@PedroA not unfamiliar, https://en.uncyclopedia.co/wiki/AAAAAAAAA!

I think it's because when the scientist heard about the virus he/she went "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa".

The question is EXACTLY how large viral bioinformatics projects are initiated, involving lots of different approaches. What is sought is a genetic observation that is anomalous and assess a correlation with virological phenotypes. In this instance, I wouldn't embark on this project, but someone else might disagree. For example, a long time ago we spotted a given untranslated region of a given virus (not Cov) was involved in attenuation (vaccine strain anomalies). Today the concept is a designed vaccine in stage III trials looking good. When bioinformatians say stuff, other virologists listen

The scientist must have died while sequencing it?

@NONONO the NOP-sled analogy would make sense at the start of the gene but this is the end

"Help, I'm trapped in a gene-sequencing lab!"

He died from the virus and feel onto the a button

Good observation! The 3' poly(A) tail is actually a very common feature of positive-strand RNA viruses, including coronaviruses and picornaviruses.
For coronaviruses in particular, we know that the poly(A) tail is required for replication, functioning in conjunction with the 3' untranslated region (UTR) as a cis-acting signal for negative strand synthesis and attachment to the ribosome during translation. Mutants lacking the poly(A) tail are severely compromised in replication. Jeannie Spagnolo and Brenda Hogue report:

The 3′ poly (A) tail plays an important, but as yet undefined role in Coronavirus genome replication. To further examine the requirement for the Coronavirus poly(A) tail, we created truncated poly(A) mutant defective interfering (DI) RNAs and observed the effects on replication. Bovine Coronavirus (BCV) and mouse hepatitis Coronavirus A59 (MHV-A59) DI RNAs with tails of 5 or 10 A residues were replicated, albeit at delayed kinetics as compared to DI RNAs with wild type tail lengths (>50 A residues). A BCV DI RNA lacking a poly(A) tail was unable to replicate; however, a MHV DI lacking a tail did replicate following multiple virus passages. Poly(A) tail extension/repair was concurrent with robust replication of the tail mutants. Binding of the host factor poly(A)- binding protein (PABP) appeared to correlate with the ability of DI RNAs to be replicated. Poly(A) tail mutants that were compromised for replication, or that were unable to replicate at all exhibited less in vitro PABP interaction. The data support the importance of the poly(A) tail in Coronavirus replication and further delineate the minimal requirements for viral genome propagation.
Spagnolo J.F., Hogue B.G. (2001) Requirement of the Poly(A) Tail in Coronavirus Genome Replication. In: Lavi E., Weiss S.R., Hingley S.T. (eds) The Nidoviruses. Advances in Experimental Medicine and Biology, vol 494. Springer, Boston, MA

Yu-Hui Peng et al. also report that the length of the poly(A) tail is regulated during infection:

Similar to eukaryotic mRNA, the positive-strand coronavirus genome of ~30 kilobases is 5’-capped and 3’-polyadenylated. It has been demonstrated that the length of the coronaviral poly(A) tail is not static but regulated during infection; however, little is known regarding the factors involved in coronaviral polyadenylation and its regulation. Here, we show that during infection, the level of coronavirus poly(A) tail lengthening depends on the initial length upon infection and that the minimum length to initiate lengthening may lie between 5 and 9 nucleotides. By mutagenesis analysis, it was found that (i) the hexamer AGUAAA and poly(A) tail are two important elements responsible for synthesis of the coronavirus poly(A) tail and may function in concert to accomplish polyadenylation and (ii) the function of the hexamer AGUAAA in coronaviral polyadenylation is position dependent. Based on these findings, we propose a process for how the coronaviral poly(A) tail is synthesized and undergoes variation. Our results provide the first genetic evidence to gain insight into coronaviral polyadenylation.
Peng Y-H, Lin C-H, Lin C-N, Lo C-Y, Tsai T-L, Wu H-Y (2016) Characterization of the Role of Hexamer AGUAAA and Poly(A) Tail in Coronavirus Polyadenylation. PLoS ONE 11(10): e0165077

This builds upon prior work by Hung-Yi Wu et al, which showed that the coronaviral 3' poly(A) tail is approximately 65 nucleotides in length in both genomic and sgmRNAs at peak viral RNA synthesis, and also observed that the precise length varied throughout infection. Most interestingly, they report:

Functional analyses of poly(A) tail length on specific viral RNA species, furthermore, revealed that translation, in vivo, of RNAs with the longer poly(A) tail was enhanced over those with the shorter poly(A). Although the mechanisms by which the tail lengths vary is unknown, experimental results together suggest that the length of the poly(A) and poly(U) tails is regulated. One potential function of regulated poly(A) tail length might be that for the coronavirus genome a longer poly(A) favors translation. The regulation of coronavirus translation by poly(A) tail length resembles that during embryonal development suggesting there may be mechanistic parallels.
Wu HY, Ke TY, Liao WY, Chang NY. Regulation of coronaviral poly(A) tail length during infection. PLoS One. 2013;8(7):e70548. Published 2013 Jul 29. doi:10.1371/journal.pone.0070548

It's also worth pointing out that poly(A) tails at the 3' end of RNA are not an unusual feature of viruses. Eukaryotic mRNA almost always contains poly(A) tails, which are added post-transcriptionally in a process known as polyadenylation. It should not therefore be surprising that positive-strand RNA viruses would have poly(A) tails as well. In eukaryotic mRNA, the central sequence motif for identifying a polyadenylation region is AAUAAA, identified way back in the 1970s, with more recent research confirming its ubiquity. Proudfoot 2011 is a nice review article on poly(A) signals in eukaryotic mRNA.

Maybe worth noting that the poly(A) tail in eukaryotic mRNAs serves as a binding region for Poly(A)-binding protein which helps take mRNA out of the nucleus.

User NONONO mentioned in comments a NOP sled, a trick used in computer viruses, used to raise the probability of a piece of code (randomly injected into the host) to be actually ran. In that context, it's always a starting point. Is it possible that this tail-of-As, as described here, is actually not a "tail" (end of data, last piece translated, etc), but a "head" instead? A "head" that is "caught" or "detected" by something and then it is where the process begins and proceeds with reading/processing/replication, or fails if the A-chain is too short and the process "loses its grasp" on it?

Knowing that these entities all share a similar 3' poly(A) tail and that it helps in replication, is it possible to essentially bind off the end with a synthetic molecule, or would that interfere with regular host operation?

What does the apostrophe in 3' mean?  Is that "Three feet", "Three Prime" ?

@JPhi1618 3-prime

@HouseCat Genes are translated to RNA, which have a 3' polyA tail. Doing so would prevent protein synthesis in your cells, and you don't want that. Ricin is a poison that blocks protein synthesis, though it does so through a different mechanism

@quetzalcoatl This is getting off topic, but yes, polyA tails are important to promote translation because the mRNA gets put head to tail in a circle. Here's a paper on this "circularization". https://pubmed.ncbi.nlm.nih.gov/9702200/

This question is quite general, so I'm going to attempt to tie it back to bioinformatics.
Background
The tree for the current coronavirus is here, showing it is closely related to bat-coronavirus and in particular SARS.
Question
The bioinformatics question for the current coronavirus is why this virus appears to be able to infect humans and transmit to human.
Genome size
Firstly, you said that 30kb was large ... this is a standard size for a coronavirus genome, albeit it is unusual in that the family Coronaviridae are the largest genomes for a single stranded RNA virus, for example flaviviruses are 10kb. Thus, all coronaviruses are all approximately 30Kb. Some coronaviruses don't infect humans (zero symptoms), some cause very mild symptoms, others are MERS and SARS with 40-60% and 10% mortality rates, respectively. So, genome size is of little bioinformatics interest in my opinion.
Polyadenylation
Polyadenylation and capping (5' methylation) enable the RNA to be trafficked and transcribed by ribosomes and the mechanism is widely used by viruses. Methylation would also prevent the innate immune response from the shredding the vRNA. Koonin and Moss (2010), interpreted a given capping mechanism as being common to the Mononegavirales - a viral Order including measles, mumps, Ebolavirus. Its a big statement, but regardless poly-A and capping are simply mimicking the host mRNA which a lot of viruses use. Poly-A and capping per se are not really interesting.
Evolution and SARS A more detailed examination the evolution of 2019-nCov and its epidemiology in relation to SARS can be found here
Conclusion
The bioinformatics question is the genome size wierd - no, its standard for a coronavirus, is the poly-A weird - no its generic amongst lots of viruses as is capping. Is the length of the poly-A excessive (33 As), it looks odd  but a human genecist/bioinformaticist needs to answer that ... so is it (potentially) linked with its epidemiology/clinical symptoms?
I don't think 33 poly-As are linked with anything bioinformatically interesting. This is because it will likely vary dramatically between genomes (not simply epidemic vs. non-epidemic strains). I don't know the mechanism for poly-adenylation, but I think slippage is a likely mutation resulting in large variations between individual genomes, particularly for poly-A - which notorious for slippage.
So ultimately could poly-As be linked with the ability of the new coronavirus to infect/transmit and could we therefore explore that bioinformatically? I personally think slippage mutations would prevent a clonal lineage emerging, i.e. that the size of the poly-As is not stable between genomes, but that would assume a given given mechanism of polyadenylation. Thus as a bioinformatics question I wouldn't pursue it, because I don't think there is sufficient biological rationale. I agree weird stuff should be questioned and that bit of the genome jumps out ... but I doubt it would go anywhere.
Slippage
The definition of a slippage mutation is here, but basically it means this genome has 33 poly-As, however another isolate from the same epidemic could say have 30 poly-As (just an example), another might have 25 poly-As and so on.
Just my 2 cents

Can anybody give me some pointers on understanding that relationship tree? It looks like the current virus split from a bat virus. And that the same ancestral virus is also ancestral to a number of other viruses in both humans and bats. Is that correct?

@puppetsock this is separate question and the rule of BioSE is one question per post, otherwise it will be moderated. If you post the question I will provide you with a detailed answer, with additional detail about the shared SARS receptor

Some of the other answers here seem quite good; at the same time I think the core answer to the OP's question is maybe a bit hard to tease out of them, so I'd like to try to state it more plainly. It's worth noting that a truly complete answer to this question seems to be beyond current research, but any kind of "Why?" is inevitably a hard or even impossible sort of question to answer fully in biology. We have some ideas about it though.
mRNA is used as a template for protein synthesis within a cell. A single mRNA is used repeatedly, but is eventually "used up" and taken apart. In eukaryotes, poly(A) tails are almost always found on mRNAs produced in the nucleus. The poly(A) tail is ultimately shortened during the transcription process, and this shortening contributes to the mRNA being degraded. (See here for more.)
Coronaviruses also have a poly(A) tail, similarly to eukaryote mRNA. The precise mechanical functions of this poly(A) tail and the means of its synthesis are objects of ongoing research, but research has shown that its presence greatly increases the degree to which Coronavirus RNA is replicated by the host cell. Research has also shown that longer tails increase replication compared to shorter tails. It's quite likely that the presence of the tail assists in recruiting the cell's protein synthesis machinery and allows the RNA to last longer within the host cell, just as it does in the cell's own mRNA.
Interestingly, the pattern in which Coronavirus poly(A) tail length is regulated during infection, in which it starts out shorter, gets longer, then gets much shorter, resembles poly(A) tail length regulation of mRNA during eukaryote embryogenesis, suggesting parallels (see the paper in the "longer tails" link for more on this as well). Longer poly(A) tail length is closely tied to greater translational efficiency in that context.
There has been some speculation in the comments as to whether or not the Coronavirus poly(A) tail resembles a NOP sled in computer programming. I think the resemblance is mostly coincidental. NOP sleds are used in exploits because a processor, encountering a NOP, moves to the next instruction without taking any other actions. A long chain of NOPs, if entered by the processor at any point within it, will lead it to the instructions at the "bottom," after the NOPs. This is advantageous to use if you can't get the processor to go exactly where you want but you know it will end up somewhere close by, because it increases your chances of having your payload executed.
It's unusual to see a lengthy NOP sled in legitimate code, to the point that people writing them usually have to disguise their function in order to avoid automatic detection. (see pg. 183 here) In contrast, a poly(A) tail is almost universally found on nuclear eukaryote mRNA (and on some mRNAs of almost all organisms to some capacity, even mitochondria). Furthermore, the functions of the poly(A) tail are complex enough that it's still an object of ongoing research decades after its initial discovery, whereas a NOP sled does one very mechanical thing. Since the environment inside a cell is so different from the environment of a processor interacting with memory, I think it's hard to make comparisons that are so granular as to deal with a specific set of machine instructions, at least in this kind of context—a processor is a very straightforward kind of machine compared to a cell.

Hi Zoe, you stated "the pattern in which Coronavirus poly(A) tail length is regulated during infection, in which it starts out shorter, gets longer, then gets much shorter..". The analogy you used was Drosophila development from egg to the maturation of the oocyte. Interesting idea. Could you cite the variation of poly-A for coronavirus please? The context isn't clear whether it is within a life-cycle within a cell, within an infection (patient) or between patients (transmission). I suspect it is an in vitro observation and represents the infection of an established cell line.

It's in this paper. Specifically: "In this study we report that the coronaviral 3’-terminal poly(A) tail length in total viral RNA, sgmRNA7, and DI RNA is relatively short (~26-45 nt) in infected cells at 0-2 hpi, increases to peak length (~65 nt) at ~6-10 hpi, and gradually decreases in size (~30-45 nt) after ~10 h of infection." They also draw the analogy to embryo development in their abstract. You're right about the context—it's in vitro using confluent human adenocarcinoma (HRT-18) cells.

Thanks its an in vitro time series over ~5 days, where the poly-A tail peaks in size quite early in the time-series, albeit the decay in poly-A size thereafter is IMO not a great analysis of their northern blot and would have been better analysed by a data scientist. The key message is clear the poly-A tail length is dynamic for Cov for the in vitro model and therefore possibly during the course of a single patient infection.

@ZoëSparks, I'm interested in your feedback regarding this comment, specifically 'Is it possible that this tail-of-As, as described here, is actually not a "tail" (end of data, last piece translated, etc), but a "head" instead?'. Let's interpret 'head' as "earlier, causality-wise", since molecular biology doesn't have any 'top' or 'bottom'.

@daveloyall Well, in this case, there is a top and bottom in the sense of where translation starts and stops on the mRNA; it begins past the 5' cap and ends before the poly(A) tail, hence the terms (see this diagram for a nice illustration). As far as the tail being "caught or detected" to begin synthesis goes, it is part of what causes mRNA to be taken up, but not the only factor—the cap plays a significant role as well, for instance (see here for more).

@daveloyall As a side note, the poly(A) tail only promotes recruitment in eukaryotes; in prokaryotes, as well as eukaryote organelles like mitochondria and chloroplasts, the presence of a poly(A) tail promotes degradation of mRNA (see here for instance).

Not an expert, but some searching on eukaryotic positive-strand RNA viruses seems to show that polyadenylation is not uncommon.  For example, Steil, et al., 2010.


text:
Pulling out a certain gene in a volcano plot
I have generated a volcano plot with a differential expression file.
Code for inputing file:
macrophage_list <- read.table("differential_expression_macrophage.csv", header = T, sep = ",")
library(tidyr)
final_df <- df %>% pivot_longer(., -c(Feature.ID, Feature.Name), names_to = c("set",".value"), names_pattern = "(.+)_(.+)")

# A tibble: 80 x 6
   Feature.ID Feature.Name set       Mean.Counts Log2.fold.change Adjusted.p.value
   <fct>      <fct>        <chr>           <dbl>            <dbl>            <dbl>
 1 a          A            Cluster.1    0.000961            0.292           1
 2 a          A            Cluster.2    0.000902            0.793           1
 3 a          A            Cluster.3    0.00181             1.46            0.758
 4 a          A            Cluster.4    0.000642            0.269           1
 5 b          B            Cluster.1    0.000320            1.95            0.910
 6 b          B            Cluster.2    0.00180             4.77            0.154
 7 b          B            Cluster.3    0                   2.19            1
 8 b          B            Cluster.4    0                   1.66            1
 9 c          C            Cluster.1    0.00128            -2.01            0.0467
10 c          C            Cluster.2    0.00632             0.352           1
# … with 70 more rows

output of head(final_tumor)
# A tibble: 6 x 6
  Feature.ID        Feature.Name set       Mean.Counts Log2.fold.change Adjusted.p.value
  <fct>             <fct>        <chr>           <dbl>            <dbl>            <dbl>
1 ENSG00000227232.5 WASH7P       Cluster.1     0                   1.50            1
2 ENSG00000227232.5 WASH7P       Cluster.2     0                   1.73            1
3 ENSG00000227232.5 WASH7P       Cluster.3     0                   1.77            1
4 ENSG00000227232.5 WASH7P       Cluster.4     0.00114             4.30            0.293
5 ENSG00000227232.5 WASH7P       Cluster.5     0                   2.15            1
6 ENSG00000227232.5 WASH7P       Cluster.6     0                   1.22            1

output of tail(final_tumor)
# A tibble: 6 x 6
  Feature.ID        Feature.Name set        Mean.Counts Log2.fold.change Adjusted.p.value
  <fct>             <fct>        <chr>            <dbl>            <dbl>            <dbl>
1 ENSG00000210196.2 MT-TP        Cluster.6       0.0699          -0.202           0.790
2 ENSG00000210196.2 MT-TP        Cluster.7       0.0801           0.0386          1
3 ENSG00000210196.2 MT-TP        Cluster.8       0.0711           0.0875          1
4 ENSG00000210196.2 MT-TP        Cluster.9       0.0152          -2.31            0.00127
5 ENSG00000210196.2 MT-TP        Cluster.10      0.0147          -2.30            0.00612
6 ENSG00000210196.2 MT-TP        Cluster.11      0.122            0.762           1

Code for generating volcano plot:
library(ggplot2)
library(ggrepel)
ggplot(final_tumor, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(final_tumor, Adjusted.p.value < 0.05),
                  aes(label = Feature.Name))

Now, I want to pull out a certain gene, Casp14, from the list and box it on the plot. How do I do that?

As you did for labeling genes with an adjusted p value below 0.05, you can subset your dataset for keeping only rows corresponding to "Casp14":
library(ggplot2)
library(ggrepel)
ggplot(final_tumor, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(final_tumor, Adjusted.p.value < 0.05),
                  aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"),
                  aes(label = Feature.Name), color = "red")

With this code, you should see now the labeling of the gene of interest (Casp14) in red on your volcano plot.

How about for Enhanced volcano? I have this code EnhancedVolcano(final_tumor, lab = as.character(final_tumor$FeatureName), x = 'Log2.fold.change', y = 'Adjusted.p.value', xlim = c(-8,8), title = 'Tumor', pCutoff = 10e-5, FCcutoff = 1.5, pointSize = 3.0, labSize = 3.0)

ggplot2 is taking forever to generate.

How long is your dataset ? Maybe it is not relevant to plot all genes. Maybe, you can plot a representative part of your dataset (e.g. 2000 randomly selected genes). Regarding EnhancedVolcano, I do not have experience with it. Sorry.

How do I go about just plotting a representation of my dataset?

you can use sample function to select random rows of your dataframe: rep_final = final_tumor[sample(1:nrow(final_tumor), size = 1000), ]

Where in ggplot(final_tumor, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+    geom_point()+   geom_text_repel(data = subset(final_tumor, Adjusted.p.value < 0.05),                    aes(label = Feature.Name))+   geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"),                    aes(label = Feature.Name), color = "red") do I add it?

Do it before ggplot to create a subset of your dataframe final_tumor and then replace final_tumor by rep_final in the ggplot command

Is there a way to make it look nicer? I uploaded the image to the question.

To reduce the number of labeled genes on your graph, you could choose a more stringent cutoff for the labeling. For example, instead of labeling genes with p values below 0.05, you can plot those that have an absolute logFC higher than 1 or 2 or 3 ... whatever. At the end, it's up to you. You have to think what is important to show with this graph and what should be the message that you sent with this graph. BTW, your question in comments is beyond the question of this post ;)

Okay, I will create a new question.

Sure. As you wish. By the way, as you are a new contributor, you should take a second to read this: https://bioinformatics.stackexchange.com/help/someone-answers


text:
probability of an entire sequence, not just one value
I am reading "Statistical methods in bioinformatics"
I came across this formula:

As explained in the book O are the observations:

My question is why the author only specified $prob(O)$ for a sequence of length=1 ?, I assume he is only talking about one observation, but generally we have a sequence $O_1 ,O_2,O_3,...$.
I assume that $Prob(sequence)= P(O_1).P(O_2)...$ but I am not sure.

The capital non-indexed O seems to refer to the whole observed sequence (O = O1, O2, ...), so the probability P(O) is of observing the whole sequence of unspecified length. Or, am I missing something?

The author writes in the quote that $O$ is to be taken for the whole set of positions $O_{i}$.
In the extracts that you provide there is no reason to believe that it is only for a sequence of length 1. The author writes "the observed sequence $O = O_1, O_2...$"
Thus as you say $Pr(sequence) = Pr(O) = Pr(O_1) * Pr(O_2)...$ according to my understanding of the quotes.