# Question
Title: How to extract metadata from NCBI's short read archive (SRA) for a few runs?
I wish to extract metadata from a list of runs on NCBI's short read archive. For instance, I'd like to extract the library name ("HS0798") from the following run info:

https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=run\_browser&run=SRR390728

I have to extract metadata for a few hundred run numbers. Is there a existing command-line tool (sratoolkit, python, perl, bash, etc) that will extract the *metadata* about the run?

I could grep it out of a webpage dump, but I notice there are XML schemas provided on the site, which to me indicates that there's some sort of XML API available somewhere.

# Answer
Entrez Direct is perfect for this sort of thing. You can do something like this:

```
$ esearch -db sra -query 'SRR390728' \
  | efetch -format runinfo \
  | cut -f1,12 -d,
Run,LibraryName
SRR292241,HS0798
SRR390728,HS0798

```

The `runinfo` table format has 47 fields including experiment accession, project accessions, bioproject and biosample identifiers and such. You may want to explore that a little bit more if you are interested in downloading metadata in addition to the LibraryName.

Also, if you prefer to parse XML output, you can append `-mode xml` to the `efetch` command to retrieve the runinfo table in XML format. `xtract` tool that comes with the Entrez Direct package can then be used to convert that into a table.
> 10 votes

# Answer
Here's a very crude script version, which should work without downloading fancy software tools. It will only work as long as the website stays in its current layout.

It extracts the particular `<td>` from the webpage, based on a nearby table link ID. and prints it out.

Note: This is brittle. It will only work for library names that fit in a single-line of the HTML, and will ignore some webserver errors. The content is also not html-decoded (so its output might show things like `&` if the field contained html-special characters).

```
#!/bin/bash

set -e -o pipefail

# pass the run as a single argument to the script
run="$1"

# extract library name field. second cell in visual order.
curl -o - "https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=run_browser&run=$run" 2>/dev/null | \
   grep 'a id="exp-accession"' -A 7 | \
   tail -n +3 | head -n 1 | \
   sed -e 's@^<td>@@' -e 's@</td>$@@'

```
> 2 votes

# Answer
Its got to be the SRA Toolkit The way I would do it (a traditional Perl digger) doesn't work on the Tracer NCBI site, at least not readily and dumping the XML is likewise unlikely to succeed.
> 1 votes

# Answer
I've been told by the support team that there's a way to pull it out in csv form by GETting the following URL, replacing `<SRRNUMBER>`, with, e.g. `SRR000001`:

```
https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?save=efetch&db=sra&rettype=runinfo&term=<SRRNUMBER>

```

So this would be:

```
# download the file attachment
wget -nv -O foo.csv  "https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?save=efetch&db=sra&rettype=runinfo&term=SRR000001"

# or just get the output on stdout
curl -v "https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?save=efetch&db=sra&rettype=runinfo&term=SRR000001"

```

As an example, the output looks like this for SRR000001. It namely contains additional info on other runs from the same SRX:

```
Run,ReleaseDate,LoadDate,spots,bases,spots_with_mates,avgLength,size_MB,AssemblyName,download_path,Experiment,LibraryName,LibraryStrategy,LibrarySelection,LibrarySource,LibraryLayout,InsertSize,InsertDev,Platform,Model,SRAStudy,BioProject,Study_Pubmed_id,ProjectID,Sample,BioSample,SampleType,TaxID,ScientificName,SampleName,g1k_pop_code,source,g1k_analysis_group,Subject_ID,Sex,Disease,Tumor,Affection_Status,Analyte_Type,Histological_Type,Body_Site,CenterName,Submission,dbgap_study_accession,Consent,RunHash,ReadHash
SRR000001,2008-04-04 14:36:13,2015-04-07 21:54:15,470985,129515135,236041,274,298,,https://sra-download.ncbi.nlm.nih.gov/traces/sra27/SRR/000000/SRR000001,SRX000007,SID2748,WGS,RANDOM,GENOMIC,PAIRED,3000,0,LS454,454 GS FLX,SRP000001,PRJNA33627,1,33627,SRS000100,SAMN00001583,simple,9606,Homo sapiens,NA18505,,,,,female,,no,,,,,454MSC,SRA000197,,public,25D46C95842BCD2F235DDD48070341C4,5FC1B56C9A5332BD21A6FCC6713C3159
SRR000013,2008-04-04 15:15:44,2012-01-19 15:14:06,478304,131410350,316908,274,307,,https://sra-download.ncbi.nlm.nih.gov/traces/sra4/SRR/000000/SRR000013,SRX000007,SID2748,WGS,RANDOM,GENOMIC,PAIRED,3000,0,LS454,454 GS FLX,SRP000001,PRJNA33627,1,33627,SRS000100,SAMN00001583,simple,9606,Homo sapiens,NA18505,,,,,female,,no,,,,,454MSC,SRA000197,,public,A576F73351E329617F91C92845ADE3BE,0EFE27254FB8790C80E4A8A3A3CCBF1D
SRR000019,2008-04-04 15:20:15,2012-01-19 15:14:06,537808,147928455,358612,275,347,,https://sra-download.ncbi.nlm.nih.gov/traces/sra4/SRR/000000/SRR000019,SRX000007,SID2748,WGS,RANDOM,GENOMIC,PAIRED,3000,0,LS454,454 GS FLX,SRP000001,PRJNA33627,1,33627,SRS000100,SAMN00001583,simple,9606,Homo sapiens,NA18505,,,,,female,,no,,,,,454MSC,SRA000197,,public,5E42E66B2972A28C3C1B714499B3CE1E,A44F729DF562E11BDA4F950F171085B3
SRR000033,2008-04-04 15:26:19,2012-01-19 15:14:16,476737,131230950,316780,275,307,,https://sra-download.ncbi.nlm.nih.gov/traces/sra4/SRR/000000/SRR000033,SRX000007,SID2748,WGS,RANDOM,GENOMIC,PAIRED,3000,0,LS454,454 GS FLX,SRP000001,PRJNA33627,1,33627,SRS000100,SAMN00001583,simple,9606,Homo sapiens,NA18505,,,,,female,,no,,,,,454MSC,SRA000197,,public,8D775F0FB1E78F2818166B4E6C7BC320,6C794108B0172AE98323C9A809BE892E
SRR000047,2008-04-04 15:34:57,2012-01-19 15:14:25,322279,88301883,214744,273,208,,https://sra-download.ncbi.nlm.nih.gov/traces/sra4/SRR/000000/SRR000047,SRX000007,SID2748,WGS,RANDOM,GENOMIC,PAIRED,3000,0,LS454,454 GS FLX,SRP000001,PRJNA33627,1,33627,SRS000100,SAMN00001583,simple,9606,Homo sapiens,NA18505,,,,,female,,no,,,,,454MSC,SRA000197,,public,B4D7BC3DAD22AE58EADC18B2F05B2732,154C9532706F955A63CF685C4D71D9E6
SRR000049,2008-04-04 15:35:58,2012-01-19 15:14:34,506408,139029601,328595,274,336,,https://sra-download.ncbi.nlm.nih.gov/traces/sra4/SRR/000000/SRR000049,SRX000007,SID2748,WGS,RANDOM,GENOMIC,PAIRED,3000,0,LS454,454 GS FLX,SRP000001,PRJNA33627,1,33627,SRS000100,SAMN00001583,simple,9606,Homo sapiens,NA18505,,,,,female,,no,,,,,454MSC,SRA000197,,public,60D02918CD7686DA084607BFD9D8E6A5,E9A732F2A8F36B8EB1CD38CA41BA8FB5
SRR000061,2008-04-04 15:41:48,2012-01-19 15:14:37,532884,146330433,351541,274,348,,https://sra-download.ncbi.nlm.nih.gov/traces/sra4/SRR/000000/SRR000061,SRX000007,SID2748,WGS,RANDOM,GENOMIC,PAIRED,3000,0,LS454,454 GS FLX,SRP000001,PRJNA33627,1,33627,SRS000100,SAMN00001583,simple,9606,Homo sapiens,NA18505,,,,,female,,no,,,,,454MSC,SRA000197,,public,8F6E1CA597BCAA19A852036775FA1FAB,2C035CCBD62808A01CEF6EB13BAEB40F

```

Providing an invalid SRR number will generate an empty HTTP 200 response.
> 1 votes

# Answer
Not a command-line tool but also useful is redesigned SRA Run Selector. Insert comma separated accessions and you can filter/download based on metadata.
> 0 votes

# Answer
You can get metadata for an entire project using pysradb:

```
> pysradb metadata SRP063732 --detailed

study_accession experiment_accession experiment_title                                                       experiment_desc                                                        organism_taxid  organism_name   library_strategy library_source library_selection sample_accession sample_title instrument           total_spots total_size   run_accession run_total_spots run_total_bases run_alias             sra_url                                                                                experiment_alias isolate                 dev_stage sex   tissue             treatment     BioSampleModel
 SRP063732       SRX1249861           Bisulfite-seq of rhesus macaque 26168_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 26168_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071080                    Illumina HiSeq 2000  89051187    11830170462  SRR2420581    89051187        17988339774     26168_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420581/SRR2420581.1                   26168_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249860           Bisulfite-seq of rhesus macaque 26148_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 26148_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071081                    Illumina HiSeq 2000  87335402    11763248173  SRR2420580    87335402        17641751204     26148_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420580/SRR2420580.1                   26148_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249859           Bisulfite-seq of rhesus macaque 25822_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 25822_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071082                    Illumina HiSeq 2000  83544130    9481342835   SRR2420579    83544130        13784781450     25822_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420579/SRR2420579.1                   25822_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249858           Bisulfite-seq of rhesus macaque 25154_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 25154_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071083                    Illumina HiSeq 2000  115476876   13003165573  SRR2420578    115476876       19053684540     25154_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420578/SRR2420578.1                   25154_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249857           Bisulfite-seq of rhesus macaque 23784_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 23784_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071084                    Illumina HiSeq 2000  101539017   13697691917  SRR2420577    101539017       20510881434     23784_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420577/SRR2420577.1                   23784_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249856           Bisulfite-seq of rhesus macaque 23779_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 23779_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071085                    Illumina HiSeq 2000  98394582    13326784419  SRR2420576    98394582        19875705564     23779_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420576/SRR2420576.1                   23779_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249855           Bisulfite-seq of rhesus macaque 23773_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 23773_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071086                    Illumina HiSeq 2000  84442355    11218164664  SRR2420575    84442355        17057355710     23773_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420575/SRR2420575.1                   23773_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249854           Bisulfite-seq of rhesus macaque 23764_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 23764_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071087                    Illumina HiSeq 2000  87197006    11291883405  SRR2420574    87197006        17613795212     23764_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420574/SRR2420574.1                   23764_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249853           Bisulfite-seq of rhesus macaque 23762_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 23762_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071088                    Illumina HiSeq 2000  89455608    11917710767  SRR2420573    89455608        18070032816     23762_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420573/SRR2420573.1                   23762_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249852           Bisulfite-seq of rhesus macaque 25811_LowDrinker: Nucleus accumbens    Bisulfite-seq of rhesus macaque 25811_LowDrinker: Nucleus accumbens    9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071089                    Illumina HiSeq 2000  88216999    10285962574  SRR2420572    88216999        14555804835     25811_NaC_LMD _files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420572/SRR2420572.1                   25811_NucleusAccumbens  Adult     male  Nucleus accumbens  LightEthanol  Model organism or animal
 SRP063732       SRX1249851           Bisulfite-seq of rhesus macaque 25790_LowDrinker: Nucleus accumbens    Bisulfite-seq of rhesus macaque 25790_LowDrinker: Nucleus accumbens    9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071090                    Illumina HiSeq 2000  94731241    10748979142  SRR2420571    94731241        15630654765     25790_NaC_LMD _files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420571/SRR2420571.1                   25790_NucleusAccumbens  Adult     male  Nucleus accumbens  LightEthanol  Model organism or animal
 SRP063732       SRX1249850           Bisulfite-seq of rhesus macaque 25742_LowDrinker: Nucleus accumbens    Bisulfite-seq of rhesus macaque 25742_LowDrinker: Nucleus accumbens    9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071091                    Illumina HiSeq 2000  113512952   17943763654  SRR2420570    113512952       22929616304     25742_NaC_LMD _files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420570/SRR2420570.1                   25742_NucleusAccumbens  Adult     male  Nucleus accumbens  LightEthanol  Model organism or animal
 SRP063732       SRX1249849           Bisulfite-seq of rhesus macaque 25240_LowDrinker: Nucleus accumbens    Bisulfite-seq of rhesus macaque 25240_LowDrinker: Nucleus accumbens    9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071092                    Illumina HiSeq 2000  81002663    9232986327   SRR2420569    81002663        13365439395     25240_NaC_LMD _files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420569/SRR2420569.1                   25240_NucleusAccumbens  Adult     male  Nucleus accumbens  LightEthanol  Model organism or animal
 SRP063732       SRX1249848           Bisulfite-seq of rhesus macaque 25184_LowDrinker: Nucleus accumbens    Bisulfite-seq of rhesus macaque 25184_LowDrinker: Nucleus accumbens    9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071093                    Illumina HiSeq 2000  110228135   12594727479  SRR2420568    110228135       18187642275     25184_NaC_LMD _files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420568/SRR2420568.1                   25184_NucleusAccumbens  Adult     male  Nucleus accumbens  LightEthanol  Model organism or animal
 SRP063732       SRX1249847           Bisulfite-seq of rhesus macaque 25157_LowDrinker: Nucleus accumbens    Bisulfite-seq of rhesus macaque 25157_LowDrinker: Nucleus accumbens    9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071094                    Illumina HiSeq 2000  97249804    11026452312  SRR2420567    97249804        16046217660     25157_NaC_LMD _files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420567/SRR2420567.1                   25157_NucleusAccumbens  Adult     male  Nucleus accumbens  LightEthanol  Model organism or animal
 SRP063732       SRX1249846           Bisulfite-seq of rhesus macaque 21119_LowDrinker: Nucleus accumbens    Bisulfite-seq of rhesus macaque 21119_LowDrinker: Nucleus accumbens    9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071095                    Illumina HiSeq 2000  68899204    9241244635   SRR2420566    68899204        13917639208     21119_NaC_LMD _files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420566/SRR2420566.1                   21119_NucleusAccumbens  Adult     male  Nucleus accumbens  LightEthanol  Model organism or animal
 SRP063732       SRX1249845           Bisulfite-seq of rhesus macaque 26104_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 26104_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071096                    Illumina HiSeq 2000  84859783    9604170375   SRR2420565    84859783        14001864195     26104_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420565/SRR2420565.1                   26104_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal
 SRP063732       SRX1249844           Bisulfite-seq of rhesus macaque 26082_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 26082_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071097                    Illumina HiSeq 2000  73346886    8248399265   SRR2420564    73346886        12102236190     26082_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420564/SRR2420564.1                   26082_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal
 SRP063732       SRX1249843           Bisulfite-seq of rhesus macaque 25207_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 25207_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071098                    Illumina HiSeq 2000  77306973    12235827725  SRR2420563    77306973        15616008546     25207_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420563/SRR2420563.1                   25207_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal
 SRP063732       SRX1249842           Bisulfite-seq of rhesus macaque 25407_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 25407_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071099                    Illumina HiSeq 2000  83985968    13355795218  SRR2420562    83985968        16965165536     25407_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420562/SRR2420562.1                   25407_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal
 SRP063732       SRX1249841           Bisulfite-seq of rhesus macaque 24818_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 24818_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071100                    Illumina HiSeq 2000  72782928    8484613123   SRR2420561    72782928        12009183120     24818_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420561/SRR2420561.1                   24818_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal
 SRP063732       SRX1249840           Bisulfite-seq of rhesus macaque 24698_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 24698_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071101                    Illumina HiSeq 2000  79414016    12627204820  SRR2420560    79414016        16041631232     24698_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420560/SRR2420560.1                   24698_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal
 SRP063732       SRX1249839           Bisulfite-seq of rhesus macaque 24340_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 24340_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071102                    Illumina HiSeq 2000  84772576    13436053136  SRR2420559    84772576        17124060352     24340_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420559/SRR2420559.1                   24340_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal
 SRP063732       SRX1249838           Bisulfite-seq of rhesus macaque 22770_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 22770_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071103                    Illumina HiSeq 2000  83310630    13224887149  SRR2420558    83310630        16828747260     22770_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420558/SRR2420558.1                   22770_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal

```
> 0 votes

---
Tags: phylogenetics, sratoolkit
---


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

# Answer
**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
```

> 3 votes

---
Tags: r, microarray
---


title: stacked bargraph with error bars

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

```

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

---
Tags: ggplot2
---


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

# Answer
> 4 votes
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.

---
Tags: r, ggplot2
---


title: Dodged bar chart

text:
# Question
Title: Dodged bar chart
I want to create a dodged barchart with different cell cycles for each cell type within each condition.

**dput(DF\_count)**

```
structure(list(orig.ident = c("treated", "treated", "treated",
"treated", "treated", "treated", "treated", "treated", "treated",
"treated", "treated", "treated", "treated", "treated", "treated",
"treated", "treated", "treated", "treated", "treated", "treated",
"treated", "treated", "treated", "treated", "treated", "treated",
"treated", "treated", "treated", "treated", "treated", "treated",
"untreated", "untreated", "untreated", "untreated", "untreated",
"untreated", "untreated", "untreated", "untreated", "untreated",
"untreated", "untreated", "untreated", "untreated", "untreated",
"untreated", "untreated", "untreated", "untreated", "untreated",
"untreated", "untreated", "untreated", "untreated", "untreated",
"untreated", "untreated", "untreated", "untreated", "untreated",
"untreated", "untreated", "untreated"), Phase = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L), .Label = c("G1", "G2M", "S"), class = "factor"), old.ident = structure(c(1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 1L, 2L, 3L, 4L, 5L,
6L, 7L, 8L, 9L, 10L, 11L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L,
10L, 11L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 1L, 2L,
3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 1L, 2L, 3L, 4L, 5L, 6L,
7L, 8L, 9L, 10L, 11L), .Label = c("Fibroblast", "T cell", "Macrophage",
"Stellate", "Acinar", "Endothelial", "Tumor", "B cell", "Mast cell",
"Ductal", "Islets of Langerhans"), class = "factor"), n = c(3257L,
666L, 560L, 446L, 91L, 464L, 192L, 506L, 162L, 207L, 90L, 603L,
608L, 238L, 98L, 22L, 110L, 215L, 95L, 63L, 97L, 10L, 2191L,
618L, 335L, 349L, 35L, 294L, 250L, 231L, 142L, 153L, 16L, 2910L,
1259L, 853L, 559L, 1134L, 348L, 306L, 133L, 116L, 12L, 192L,
706L, 1186L, 374L, 112L, 616L, 106L, 277L, 32L, 47L, 6L, 27L,
1715L, 1312L, 575L, 390L, 1036L, 250L, 252L, 91L, 66L, 8L, 68L
), Freq = c(53.8258139150554, 35.2008456659619, 49.4263018534863,
49.9440089585666, 61.4864864864865, 53.4562211981567, 29.2237442922374,
60.8173076923077, 44.141689373297, 45.2954048140044, 77.5862068965517,
9.96529499256321, 32.1353065539112, 21.0061782877317, 10.9742441209406,
14.8648648648649, 12.6728110599078, 32.7245053272451, 11.4182692307692,
17.1662125340599, 21.2253829321663, 8.62068965517241, 36.2088910923814,
32.6638477801269, 29.567519858782, 39.0817469204927, 23.6486486486486,
33.8709677419355, 38.0517503805175, 27.7644230769231, 38.6920980926431,
33.4792122538293, 13.7931034482759, 54.5863815419246, 33.5107798775619,
47.3362930077692, 52.6861451460886, 40.7035175879397, 49.4318181818182,
36.6467065868263, 51.953125, 50.6550218340611, 46.1538461538462,
66.8989547038328, 13.2432939410992, 31.5677402182592, 20.7547169811321,
10.5560791705938, 22.1105527638191, 15.0568181818182, 33.1736526946108,
12.5, 20.5240174672489, 23.0769230769231, 9.40766550522648, 32.1703245169762,
34.9214799041789, 31.9089900110988, 36.7577756833176, 37.1859296482412,
35.5113636363636, 30.1796407185629, 35.546875, 28.82096069869,
30.7692307692308, 23.6933797909408)), class = c("grouped_df",
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -66L), groups = structure(list(
    orig.ident = c("treated", "treated", "treated", "treated",
    "treated", "treated", "treated", "treated", "treated", "treated",
    "treated", "untreated", "untreated", "untreated", "untreated",
    "untreated", "untreated", "untreated", "untreated", "untreated",
    "untreated", "untreated"), old.ident = structure(c(1L, 2L,
    3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 1L, 2L, 3L, 4L, 5L,
    6L, 7L, 8L, 9L, 10L, 11L), .Label = c("Fibroblast", "T cell",
    "Macrophage", "Stellate", "Acinar", "Endothelial", "Tumor",
    "B cell", "Mast cell", "Ductal", "Islets of Langerhans"), class = "factor"),
    .rows = list(c(1L, 12L, 23L), c(2L, 13L, 24L), c(3L, 14L,
    25L), c(4L, 15L, 26L), c(5L, 16L, 27L), c(6L, 17L, 28L),
        c(7L, 18L, 29L), c(8L, 19L, 30L), c(9L, 20L, 31L), c(10L,
        21L, 32L), c(11L, 22L, 33L), c(34L, 45L, 56L), c(35L,
        46L, 57L), c(36L, 47L, 58L), c(37L, 48L, 59L), c(38L,
        49L, 60L), c(39L, 50L, 61L), c(40L, 51L, 62L), c(41L,
        52L, 63L), c(42L, 53L, 64L), c(43L, 54L, 65L), c(44L,
        55L, 66L))), row.names = c(NA, -22L), class = c("tbl_df",
"tbl", "data.frame"), .drop = TRUE))

```

Used this code to generate a stacked bar graph:

```
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)

```

But want the 2 conditions to be next to each other on the same graph. How do I use this dataset and create a dodged bar chart to have both conditions together?

# Answer
Here is a rather straight-forward work-around for your problem:

```
ggplot(DF_Count, aes(x = orig.ident, y = Freq, fill = Phase))+
  geom_col()+
  geom_text(aes(label = paste(round(Freq, 2),"%")),position = position_stack(vjust = 0.5))+
  facet_wrap(~old.ident)

```

It does not give a single plot or a dodged bar chart representing the two conditions you would like to visualize. For an exact solution to your problem, see this (actually the answer below the linked answer is analogous to what I have suggested above).
> 1 votes

---
Tags: r, ggplot2
---


title: Bar plot with error bars

text:
# Question
Title: Bar plot with error bars
I am trying to plot a box plot/stacked bar graph with error bars.

**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))`

```

Gives me the following dataset.

```
       Var1      Var2 Freq
1   Fibroblast   treated 6051
2       T cell   treated 1892
3   Macrophage   treated 1133
4     Stellate   treated  893
5       Acinar   treated  148
6  Endothelial   treated  868
7   Fibroblast untreated 5331
8       T cell untreated 3757
9   Macrophage untreated 1802
10    Stellate untreated 1061
11      Acinar untreated 2786
12 Endothelial untreated  704

library(dplyr)
DF <- df %>% rename(Cell_Type = Var1, Condition = Var2) %>%
  group_by(Condition) %>%
  mutate(Percent = Freq / sum(Freq)*100)

# A tibble: 12 x 4
# Groups:   Condition [2]
   Cell_Type   Condition  Freq Percent
   <fct>       <fct>     <int>   <dbl>
 1 Fibroblast  treated    6051   55.1
 2 T cell      treated    1892   17.2
 3 Macrophage  treated    1133   10.3
 4 Stellate    treated     893    8.13
 5 Acinar      treated     148    1.35
 6 Endothelial treated     868    7.90
 7 Fibroblast  untreated  5331   34.5
 8 T cell      untreated  3757   24.3
 9 Macrophage  untreated  1802   11.7
10 Stellate    untreated  1061    6.87
11 Acinar      untreated  2786   18.0
12 Endothelial untreated   704    4.56

```

Calculated SD/Mean with this:

```
my_sum_1 <- DF %>%
    group_by(Cell_Type) %>%
    summarise(
        n=n(),
        mean=mean(Freq),
        sd=sd(Freq)
    ) %>%
    mutate( se=sd/sqrt(n))  %>%
    mutate( ic=se * qt((1-0.05)/2 + .5, n-1))

```

Then plotted with this:

```
ggplot(my_sum_1) +
    geom_bar( aes(x=Cell_Type, y=mean), stat="identity", fill="forestgreen", alpha=0.5) +
    geom_errorbar( aes(x=Cell_Type, ymin=mean-ic, ymax=mean+ic), width=0.4, colour="orange", alpha=0.9, size=1.5) +
    ggtitle("Cell Types")

```

I want to plot based on condition. Not sure how I would do that. I tried the `fill` command but doesn't seem to work as the object is not found, I can't seem to add it back to the data set. How do I include the condition in the bar chart?

When I'm calculating sd, I think it's combining the cell type value into one, rather than taking it individually. How do I change that?

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

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

```

# Answer
> 0 votes
Your issue is coming by the fact that when you calculate `mean`, `se`, ... you removed the `Condition` column, here is your `my_sum_1` results:

```
# A tibble: 7 x 6
  Cell_Type                n  mean    sd    se     ic
  <fct>                <int> <dbl> <dbl> <dbl>  <dbl>
1 Fibroblast               2 5691   509. 360    4574.
2 T cell                   2 2824. 1319. 932.  11849.
3 Macrophage               2 1468.  473. 334.   4250.
4 Tumor                    2  746   126.  89.0  1131.
5 Islets of Langerhans     2  202.  121.  85.5  1086.
6 Endothelial              2  786   116.  82    1042.
7 B cell                   2  544   407. 288    3659.

```

To calculate the `mean`, `sd`, ... from each cell types in each conditions, you need to group by `Cell_Type` and `Condition` by writing:

```
my_sum_1 <- DF %>%
  group_by(Cell_Type, Condition) %>%
  summarise(
    n=n(),
    mean=mean(Freq),
    sd=sd(Freq)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1))

```

But now, your table looks like:

```
# A tibble: 14 x 7
# Groups:   Cell_Type [7]
   Cell_Type            Condition     n  mean    sd    se    ic
   <fct>                <fct>     <int> <dbl> <dbl> <dbl> <dbl>
 1 Fibroblast           treated       1  6051    NA    NA    NA
 2 Fibroblast           untreated     1  5331    NA    NA    NA
 3 T cell               treated       1  1892    NA    NA    NA
 4 T cell               untreated     1  3757    NA    NA    NA
 5 Macrophage           treated       1  1133    NA    NA    NA
 6 Macrophage           untreated     1  1802    NA    NA    NA
 7 Tumor                treated       1   657    NA    NA    NA
 8 Tumor                untreated     1   835    NA    NA    NA
 9 Islets of Langerhans treated       1   116    NA    NA    NA
10 Islets of Langerhans untreated     1   287    NA    NA    NA
11 Endothelial          treated       1   868    NA    NA    NA
12 Endothelial          untreated     1   704    NA    NA    NA
13 B cell               treated       1   832    NA    NA    NA
14 B cell               untreated     1   256    NA    NA    NA

```

Because in your initial dataset, you only have one value for each cell type in each condition, so you can calculate `sd`, `sem`. You need to first merge your replicates into your dataframe in order to be able to calculate your `sem`.

Then, you will be able to do:

```
library(ggplot2)
ggplot(my_sum_1, aes(x=Cell_Type, y=mean, fill = Condition)) +
  geom_col(alpha=0.5) +
  geom_errorbar( aes(ymin=mean-ic, ymax=mean+ic), width=0.4, colour="orange", alpha=0.9, size=1.5) +
  ggtitle("Cell Types")

```

Hope it helps you to figure it out what is your issue with your current data.

---
Tags: r, ggplot2
---


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

# Answer
> 163 votes
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.

# Answer
> 32 votes
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

# Answer
> 21 votes
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.

# Answer
> 13 votes
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*.

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


# Question
Title: Visualizing the difference among three groups
I have 3 data frames for three groups of patients 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   total
1:    LP6008337-DNA_H06  927  773 40756   42456
2:    LP6008334-DNA_D02 1049  799 31009   32857
> dim(dat1)
[1] 21  6
>

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

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

> dim(dat3)
[1] 106   6

```

I want to show if the number of each mutational category and total is different between these three groups by chi-square test

Actual question is if total number of mutations, SNP, DEL and INS are statistically significant among groups. I used pairwise t test but I afraid this test is not a suitable test for the distribution of data moreover I don't know how to visualize the p-value

Imagining the comparison of total number of mutations between three groups, this picture is a good example

Or in this plot they compare several features but only two groups

Can you help ?

Thanks

# Answer
Here a way to do it is to start first by creating a dataset containing your three different group, based on the few lines you display:

```
dat1$Condition = "dat1"
dat2$Condition = "dat2"
dat3$Condition = "dat3"
colnames(dat1)[1] = "Patient"
DF <- rbind(dat1,dat2,dat3)

             Patient   DEL   INS    SNP  total Condition
1: LP6008337-DNA_H06   927   773  40756  42456      dat1
2: LP6008334-DNA_D02  1049   799  31009  32857      dat1
3: LP6008031-DNA_E01 13552  3374  62105  79031      dat2
4: LP6005500-DNA_G01   539   500  43451  44490      dat2
5: LP6005935-DNA_F03 39168 16739  58095 114002      dat3
6: LP6008269-DNA_D08   849   910 103501 105271      dat3

```

Then, you can reshape `DF` in order to put all values DEL, INS, SNP and Total into a single column and create a categorical column. Doing this will be helpful for later producing the facets of your plot. To do that, you can use `pivot_longer` from `tidyr` package:

```
library(tidyr)
DF <- DF %>% pivot_longer(.,cols = c(DEL,INS,SNP, total), names_to = "var", values_to = "val")

# A tibble: 24 x 4
   Patient           Condition var     val
   <chr>             <chr>     <chr> <int>
 1 LP6008337-DNA_H06 dat1      DEL     927
 2 LP6008337-DNA_H06 dat1      INS     773
 3 LP6008337-DNA_H06 dat1      SNP   40756
 4 LP6008337-DNA_H06 dat1      total 42456
 5 LP6008334-DNA_D02 dat1      DEL    1049
 6 LP6008334-DNA_D02 dat1      INS     799
 7 LP6008334-DNA_D02 dat1      SNP   31009
 8 LP6008334-DNA_D02 dat1      total 32857
 9 LP6008031-DNA_E01 dat2      DEL   13552
10 LP6008031-DNA_E01 dat2      INS    3374
# … with 14 more rows

```

Now, you can plot your data and use the function `stat_compare_means` from `ggpubr` package to get the significance display on each tab. Based on the size of your group, you may want to adjust for the type of statistical test to be used (take a look at `?stat_compare_means` to find the statistical test adequate to your dataset). You can pass various conditions to test as a list by using `comparisons` argument.

For facetting your graph, you can use `facet_grid` (or`facet_wrap`, try the one you preferred). Altogether, you can get the plot by writing:

```
library(ggplot2)
library(ggpubr)
ggplot(DF, aes(x = Condition, y = val, fill = Condition))+
  geom_boxplot()+
  stat_compare_means(comparisons = list(c("dat1","dat2"), c("dat1","dat3"), c("dat2","dat3")))+
  stat_compare_means()+
  facet_grid(var~., scales = "free")

```

Hope it helps you to figure it out how to represent your data.
> 3 votes

---
Tags: r, ggplot2
---


# Question
Title: how to calculate coverage of each base of mapped reads?
I have a RNA aligned file without any header info. I need to calculate the coverage of each mapped base for per chromosome. Is it possible to do it with perl ?

```
NS500377:28:H16A7BGXX:2:21205:16038:7658        99      chr2   3000482 255     46M     =       3000502 56      GTCCTACAGTGTGCATTTCTCGTTTTCCACGTTTCT    AAAA)FFFFFFFFFAFFFFFFFFFFF.F.FFFFFAF    XA:i:2  MD:Z:21T4T7T1   NM:i:3
NS500377:28:H16A7BGXX:3:11503:24953:10427       99      chr2   3000482 255     50M     =       3000513 71      GTCCTACAGTGTGCATTTCTCATTTTTCAAGTTTTTCAGT        AAAAAFFFFFFFFFFFFAFFFFFF7AAFFFFFFFFFFFFF
        XA:i:1  MD:Z:21T7C10    NM:i:2
NS500377:28:H16A7BGXX:3:12404:23929:9392        163     chr2   3000482 255     49M     =       3000483 40      GTCCTACAGTGTGCATTTCCCTTTTTTCACGTTTTTTAG AAFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF XA:i:1  MD:Z:19T16C2    NM:i:2
NS500377:28:H16A7BGXX:3:12404:23929:9392        83      chr2   3000483 255     49M     =       3000482 -40     TCCTACAGTGTGCATTTCCCTTTTTTCACGTTTTTTAGT AFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF XA:i:2  MD:Z:18T16C3    NM:i:2
NS500377:28:H16A7BGXX:3:21606:15675:1639        83      chr2   3000483 255     50M     =       3000480 -43     TCTTACAGTGTGCATTTCTCTTTTTTCACGTTTTTCAGTG        AFFFFFFFFFAFFFFFFFFFFFFFFFFFFFFFFFFAAAAA
        XA:i:0  MD:Z:2C37       NM:i:1
NS500377:28:H16A7BGXX:4:12605:18505:11088       83      chr2   3000486 255     50M     =       3000481 -45     TACAATGTGCATTTCTCTTTTTTCACGTTTTTCAGTATTT        FFFFFFFAFFF)FFFFFFFFF<FFFF<AFFFFFFFAA<AA
        XA:i:1  MD:Z:4G31G3     NM:i:2
NS500377:28:H16A7BGXX:2:13303:20634:12085       99      chr2   3000487 255     49M     =       3000521 73      ACAGTGTACATTTCTCATTTTTCACGTTTTTCAGTGTTA AAAAAFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF. XA:i:2  MD:Z:7G8T21T0   NM:i:3

```

# Answer
> 1 votes
You are looking for `bedtools genomecov`, specifically `-d` option. See https://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html , for example:

```
bedtools genomecov -i A.bed -g my.genome -d

```

---

I recommend not reinventing the wheel: use `bedtools` instead of reimplementing the same functionality in Perl.

---
Tags: rna, coverage, awk
---


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

# Answer
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`:

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

---
Tags: fasta, ngs, variant-calling, gatk, qc
---


title: Kallisto for Nanopore?

text:
# Question
Title: Kallisto for Nanopore?
Although the Kallisto software was designed for RNA, it has also been very reliable in counting selected FASTA sequences in DNA sequencing. For example, given a WGS sample and my FASTA file. Kallisto is able to accurately return the counts for each sequence in my FASTA file.

Now, I'd like to do something similar on Nanopore long reads. Is that a similar software?

# Answer
It seems NanoCount can do what you are looking for.:

> EM based transcript abundance from nanopore reads mapped to a transcriptome with minimap2 Python package adapted from https://github.com/jts/nanopore-rna-analysis by Jared Simpson
>
> NanoCount estimates transcript abundance from ONT direct-RNA Sequencing reads mapped to a transcriptome. It uses an expectation-maximization approach like RSEM, Kallisto, salmon, etc to handle multi-mapping reads. The reads must be mapped to the transcript set using minimap2.
> 1 votes

---
Tags: nanopore
---


# Question
Title: Interpretion of my coronavirus 2019-nCov, Wuhan, China BLAST tree?
This is the BLAST tree of the latest coronavirus out of China (from Wuhan Institute of Virology, China). It seems strange that there is so much divergence from all the other coronaviruses. Is this expected of new diseases?

# Answer
> 4 votes
The evolutionary related group (clade) of betacoronaviruses you have identified share an amino acid homology of 85% and include SARS. I know this from the underlying tree published on BioRxiv of a broader group of betacoronaviruses, i.e. your data is a defined subset of the betacoronaviruses which all share a unique, single common ancester.

Lets call this group the SARS clade.

You have performed a nucleotide blast and asked NCBI to produce a NJ tree using 2019-nCov as the reference. I can tell this by the colour coding and the scale bar at the bottom right hand corner shows the genetic distance involved is reasonably larger than 15% divergence by amino acid data. The scale bar represents the number of mutations per nucleotide.

In summary, your tree is essentially a subset of the broader genetic diversity of the one published, but there is a rooting issue.

In your tree, the majority of sequences are from the 2002 SARS epidemic and the virtually zero genetic distance between them is simply because it is a rapidly transmitted outbreak. *I didn't realize that SARS has two independent origins both initially from bats. This is quite scary.*

* 2019-nCov is an outgroup within the 'SARS clade' hence it appears at the other side of the tree, i.e. they share a more distant common ancestor
* However 2019-nCov is not the most distant common ancester within this subset of the betacoronaviruses, i.e. SARS clade, this belongs to the two bat viruses ZC45 and ZX21. The program has made a likely rooting error (below).
* Again the blast omitted the majority of the betacoronaviruses, for example MERS

**Rooting issue** the reason I suggest there is a rooting error is because the BioRxiv tree using a broad sample of the betacoronaviruses placed the bat strains ZC45 and ZX21 as outgroups to the SARS clade and 2019-nCOV was immediately within that, so in this definition 2019-nCOV is an 'ingroup' within the SARS clade, whereas your tree it is an 'outgroup'. It not a huge issue, but the location of the "root" is determined by the common ancester above it (sequences with \<85% homology) and in your tree those have been omitted.

Generally, I like the analysis otherwise and gives an insight into SARS that I hadn't previosly been aware of.

---
Tags: sequence-alignment, phylogenetics, phylogeny, covid-19, sars-cov-2
---


# Question
Title: How to obtain file offsets corresponding to range query in block-gzipped file
How can we retrieve the file offsets of the zipped-chunks in a block-gzipped file that contain records for a region/window of interest (e.g. Chr01:100-Chr02:145)?

It seems like the CLI of `samtools`, or `tabix` should support that, but I can't find a way to just get the offsets in the datafile, without having the accompanying datafile at hand. I couldn't find a formal specification for those index files either.

Perhaps one can recommend a library that can parse those index files, in python, C, or golang? It's hard to accept that there isn't a reusable block-gzipped index parser somewhere, other than the one buried deep in samtools/htslib/gatk code.

**edit**

I'd initially given a small motivating example making use of file offsets to access a subportion of file on a remote server. Commenters were quick to point out that samtools supported a series of protocols for doing exactly that: ftp, http, and s3 (which requires https, which requires that samtools be compiled with libcurl support) Thanks! That example was detracting, so I removed it.

`samtools view` does fetch a subportion of a file (potentially from a remote location), but it does a few extra things, like decompressing the chunks retrieved. With this question, I'm interested in manipulating the gz chunks directly.

# Answer
> 2 votes
So if the purpose of obtaining file offsets is to access subportions of remote files via protocols that samtools supports, as mentioned in this comment, then you can make your life easier and use samtools/tabix directly on the cli, prefixing your datafile with the corresponding scheme.

One side effect of using the CLI however, is that the blocks retrieved will be decompressed (e.g. samtools view). You'd have to re-compress them if what you're trying to do is create a smaller "bam" file.

Keep in mind that if you retrieve all zipped chunks which cover your range of interest, then you might be extracting more data than you asked for. The first and last chunk may contain records outside your region, just because their neighboring records happened to be compressed in the same chunk. You'd have to be careful to filter those downstream.

**Specification for the tabix format**

If you absolutely need a custom solution, the tabix data structure format, authored by Heng Li, is on the samtools website:

https://samtools.github.io/hts-specs/tabix.pdf

A more in depth discussion of how the data structure works can be found here:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3042176/

The tree-shaped binning index is a bit hard to understand from the text alone, but the tabix paper mentions that their approach is almost identical to the binning algorithm in this paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC186604/pdf/X6.pdf (page 1003). That paper has a nice illustration which makes things clearer.

**Operation**

Note: This is from my understanding of how the indices interoperate.

You don't need the data file to retrieve chunk offsets, you can do the work with just the tabix file. But it's complicated, and probably error prone if you're implementing this yourself.

Let's assume the query interval is `Chr01:100-Chr02:2000`.

Each sequence (Chr01, Chr02, etc) has its own pair of indices (binning index + linear index) in the tabix file. So you would break the query into two queries: `Chr01:100-999999999999` and `Chr02:0-2000`, and concatenate the two results.

Within one sequence, to find all the blocks for open interval `[X,Y)`, you would use the binning index to retrieve all the bins overlapping this interval (There is C code to do this in the data structure pdf). Then, you need to go through all the possible bins selected, and filter chunks. Not all chunks in a given bin will be overlapping with `X, Y`.

To filter chunks within a bin, you need the help of the linear index. There is one entry in the linear index per 16kbp of sequence positions. Each entry of the linear index contains the virtual offset of the leftmost record overlapping that 16kb window. I don't know what the entry contains if there is no record overlapping that 16kb window, presumably 0.

You would find the linear entry in which position X falls (X/16Kb), take that offset, call it O\_x. You find the linear entry corresponding to position Y, and take that offset, call it O\_y. When you go through the selected bins, filter out any chunk whose offset is \< O\_x, or whose offset \>= O\_y.

Assuming your datafile is sorted, the chunk range you're interested in is contained between the minimum offset and the maximum offset selected.

The virtual offsets are 64bit unsigned integers. The higher 48 bits of the virtual chunk offset corresponds to the file offset of the datafile (the other 16 you don't need. They are the offset of the record within that decompressed chunk).

**Software**

The software I have found typically mimics closely the hts\_lib API. And hts\_lib doesn't export a function which loads an index without also requiring a datafile. As far as I know.

---
Tags: sam, samtools, indexing, tabix, bgzip
---