# Question
Title: Is there a standard tool used to convert a VCF to a BEDPE?
Many popular SV callers output a VCF. Unfortunately, there isn't a unified system at the present to label events with the same notation.

However, is there a standard method for converting these VCFs to BEDPEs?

svtools comes to mind: https://github.com/hall-lab/svtools

There are also some Bioconductor packages.

However, is there currently a standard widely used? Normally this is because a certain consortium did it in a certain way.

EDIT: There was a question in the comments about VCF vs BEDPE for structural variants.

https://github.com/jaquol/blog-like-old/blob/master/01-How-interpret-structural-variants-BEDPE-format.md

This is a structural variant in VCF format:

```
##fileformat=VCFv4.2
##source=LUMPY
##INFO=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variant">
##INFO=<ID=SVLEN,Number=.,Type=Integer,Description="Difference in length between REF and ALT alleles">
##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">
##INFO=<ID=STRANDS,Number=.,Type=String,Description="Strand orientation of the adjacency in BEDPE format (DEL:+-, DUP:-+, INV:++/--)">
##INFO=<ID=IMPRECISE,Number=0,Type=Flag,Description="Imprecise structural variation">
##INFO=<ID=CIPOS,Number=2,Type=Integer,Description="Confidence interval around POS for imprecise variants">
##INFO=<ID=CIEND,Number=2,Type=Integer,Description="Confidence interval around END for imprecise variants">
##INFO=<ID=CIPOS95,Number=2,Type=Integer,Description="Confidence interval (95%) around POS for imprecise variants">
##INFO=<ID=CIEND95,Number=2,Type=Integer,Description="Confidence interval (95%) around END for imprecise variants">
##INFO=<ID=MATEID,Number=.,Type=String,Description="ID of mate breakends">
##INFO=<ID=EVENT,Number=1,Type=String,Description="ID of event associated to breakend">
##INFO=<ID=SECONDARY,Number=0,Type=Flag,Description="Secondary breakend in a multi-line variants">
##INFO=<ID=SU,Number=.,Type=Integer,Description="Number of pieces of evidence supporting the variant across all samples">
##INFO=<ID=PE,Number=.,Type=Integer,Description="Number of paired-end reads supporting the variant across all samples">
##INFO=<ID=SR,Number=.,Type=Integer,Description="Number of split reads supporting the variant across all samples">
##INFO=<ID=EV,Number=.,Type=String,Description="Type of LUMPY evidence contributing to the variant call">
##INFO=<ID=PRPOS,Number=.,Type=String,Description="LUMPY probability curve of the POS breakend">
##INFO=<ID=PREND,Number=.,Type=String,Description="LUMPY probability curve of the END breakend">
##ALT=<ID=DEL,Description="Deletion">
##ALT=<ID=DUP,Description="Duplication">
##ALT=<ID=INV,Description="Inversion">
##ALT=<ID=DUP:TANDEM,Description="Tandem duplication">
##ALT=<ID=INS,Description="Insertion of novel sequence">
##ALT=<ID=CNV,Description="Copy number variable region">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=SU,Number=1,Type=Integer,Description="Number of pieces of evidence supporting the variant">
##FORMAT=<ID=PE,Number=1,Type=Integer,Description="Number of paired-end reads supporting the variant">
##FORMAT=<ID=SR,Number=1,Type=Integer,Description="Number of split reads supporting the variant">
##FORMAT=<ID=BD,Number=1,Type=Integer,Description="Amount of BED evidence supporting the variant">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  3128_T47D
chr1    1584528 11  N   <DUP>   .   .   SVTYPE=DUP;STRANDS=-+:6;SVLEN=66541;END=1651069;CIPOS=-449,10;CIEND=-10,503;CIPOS95=-374,0;CIEND95=0,366;IMPRECISE;SU=6;PE=6;SR=0;PRPOS=2.44491e-07,3.26458e-07,4.33697e-07,5.7608e-07,7.57136e-07,9.87322e-07,1.28796e-06,1.66553e-06,2.14862e-06,2.75607e-06,3.52645e-06,4.48654e-06,5.68684e-06,7.14253e-06,8.95608e-06,1.11311e-05,1.38146e-05,1.70081e-05,2.09074e-05,2.54614e-05,3.0947e-05,3.7325e-05,4.47552e-05,5.33515e-05,6.32491e-05,7.45332e-05,8.7438e-05,0.00010173,0.000117821,0.000135629,0.000155258,0.000177125,0.000200949,0.000227194,0.000255827,0.000286183,0.000319177,0.000353501,0.000390203,0.000428746,0.000468895,0.000511378,0.000555476,0.000600756,0.000648001,0.000695831,0.000745204,0.000795631,0.000845517,0.000896997,0.00094752,0.000998877,0.00104988,0.00110053,0.00115056,0.00120038,0.00124937,0.00129849,0.00134617,0.00139206,0.00143762,0.0014814,0.00152386,0.00156521,0.00160519,0.00164404,0.00168131,0.00171717,0.00175139,0.00178475,0.00181657,0.00184663,0.00187563,0.00190327,0.00192975,0.00195479,0.00197842,0.00200183,0.00202369,0.00204549,0.00206583,0.00208582,0.00210445,0.00212224,0.00213857,0.00215452,0.00216896,0.00218263,0.00219584,0.00220845,0.00221975,0.00223061,0.00224049,0.00224956,0.00225841,0.00226617,0.00227317,0.00227931,0.00228505,0.00228996,0.00229495,0.00229948,0.00230347,0.00230755,0.00231109,0.0023146,0.00231767,0.00232087,0.0023235,0.00232593,0.00232843,0.00233037,0.00233266,0.00233471,0.0023365,0.00233822,0.00233974,0.00234121,0.00234261,0.00234391,0.00234524,0.00234656,0.00234787,0.00234901,0.00235006,0.0023513,0.0023522,0.00235324,0.00235429,0.00235519,0.00235614,0.00235695,0.00235786,0.00235864,0.00235957,0.00236033,0.00236117,0.00236198,0.00236288,0.00236365,0.00236441,0.00236514,0.00236581,0.00236671,0.0023675,0.00236851,0.0023693,0.00237023,0.00237106,0.00237196,0.00237283,0.00237371,0.00237471,0.00237556,0.00237647,0.00237733,0.00237819,0.00237897,0.00237973,0.00238056,0.00238145,0.00238222,0.00238294,0.00238394,0.0023848,0.00238545,0.00238618,0.0023871,0.00238796,0.0023889,0.00238976,0.00239053,0.00239142,0.00239222,0.0023932,0.00239412,0.00239504,0.00239599,0.00239704,0.00239815,0.00239902,0.00239984,0.00240073,0.00240184,0.00240288,0.00240385,0.00240489,0.00240597,0.00240695,0.00240787,0.00240872,0.00240953,0.00241046,0.00241129,0.00241228,0.00241313,0.00241411,0.00241502,0.00241603,0.00241712,0.00241811,0.00241934,0.00242034,0.00242141,0.00242256,0.00242369,0.00242488,0.00242598,0.00242723,0.00242845,0.00242958,0.00243082,0.00243182,0.00243282,0.00243389,0.00243485,0.00243593,0.00243699,0.00243818,0.00243921,0.00244038,0.00244162,0.00244278,0.00244386,0.002445,0.00244621,0.00244734,0.00244858,0.00244991,0.0024512,0.00245246,0.00245379,0.00245524,0.00245665,0.00245774,0.00245913,0.00246048,0.00246179,0.00246295,0.00246419,0.00246564,0.00246697,0.00246825,0.00246954,0.0024709,0.00247228,0.00247353,0.002475,0.00247639,0.00247812,0.0024795,0.00248075,0.00248215,0.00248359,0.00248517,0.00248655,0.00248793,0.00248927,0.00249072,0.00249223,0.00249351,0.00249512,0.00249648,0.00249801,0.00249976,0.00250133,0.00250293,0.00250464,0.00250618,0.00250758,0.00250932,0.00251113,0.00251256,0.00251406,0.00251567,0.00251749,0.00251919,0.00252051,0.00252197,0.00252351,0.00252518,0.00252675,0.00252834,0.00252988,0.00253136,0.00253307,0.00253469,0.00253652,0.0025383,0.0025403,0.00254212,0.00254387,0.00254575,0.00254754,0.00254942,0.0025514,0.00255316,0.00255534,0.0025573,0.00255935,0.00256123,0.00256313,0.00256486,0.00256668,0.00256838,0.00257033,0.00257235,0.00257408,0.00257612,0.00257812,0.00258004,0.00258208,0.00258394,0.00258613,0.00258805,0.00259035,0.00259238,0.00259449,0.00259664,0.00259868,0.00260091,0.0026033,0.00260538,0.00260742,0.00260951,0.00261187,0.00261382,0.00261579,0.00261808,0.00262029,0.00262253,0.00262481,0.00262695,0.00262929,0.00263173,0.00263396,0.0026362,0.00263851,0.00264069,0.00264344,0.0026458,0.00264821,0.00265024,0.00265247,0.00265488,0.00265735,0.00265965,0.00266185,0.00266402,0.00266663,0.00266919,0.00267144,0.00267383,0.00267621,0.002679,0.00268127,0.00268393,0.00268632,0.00268891,0.00269145,0.00269389,0.00269661,0.00269912,0.00270165,0.00270404,0.0027067,0.00270939,0.00271195,0.00271455,0.00271722,0.00271969,0.00272252,0.0027254,0.002728,0.00273082,0.00273369,0.00273648,0.00273952,0.00274227,0.00274511,0.00274816,0.00275103,0.00275394,0.00275687,0.00275977,0.00276241,0.00276495,0.00276767,0.00277064,0.00277351,0.00277641,0.00277909,0.00278202,0.00278452,0.00278729,0.00278998,0.00279275,0.0027958,0.00279878,0.0028016,0.00280472,0.00280771,0.0028109,0.00281415,0.00281717,0.00282038,0.00282346,0.00282626,0.00282919,0.00283207,0.00283518,0.00283843,0.00284122,0.0028439,0.00284708,0.00284994,0.00285294,0.00285587,0.00285848,0.00286126,0.0028643,0.00286701,0.00287003,0.00287284,0.00287589,0.00287886,0.00288176,0.00288469,0.00288762,0.0028909,0.00289392,0.00289664,0.00289945,0.00290244,0.00290543,0.00290824,0.00291089,0.00291354,0.00291634,0.00291902,0.00292172,0.00292413,0.0029271,0.00292984,0.002933,0.0029355,0.00293802,0.0029405,0.00294333,0.00294613,0.00117372,0.000467584,0.000186264,7.42093e-05,2.95642e-05,1.17781e-05,4.69226e-06,1.86919e-06,7.44696e-07,2.96676e-07;PREND=3.02026e-15,4.78471e-14,7.57973e-13,1.20067e-11,1.90178e-10,3.01282e-09,4.77324e-08,7.56151e-07,1.19779e-05,0.000189743,0.00300582,0.0030039,0.00300127,0.00299837,0.00299541,0.00299255,0.00299047,0.00298793,0.00298555,0.00298267,0.00298033,0.0029782,0.00297579,0.00297295,0.00297017,0.00296769,0.00296562,0.0029632,0.00296072,0.00295816,0.00295557,0.00295308,0.00295062,0.0029479,0.0029447,0.00294211,0.00293963,0.00293734,0.00293449,0.00293173,0.00292904,0.0029266,0.00292425,0.00292122,0.00291848,0.00291567,0.00291333,0.00291077,0.00290832,0.00290537,0.00290246,0.00289957,0.00289704,0.00289385,0.00289114,0.0028884,0.00288549,0.00288284,0.00287995,0.00287726,0.00287463,0.00287144,0.00286885,0.00286569,0.0028629,0.00286073,0.00285801,0.00285472,0.00285196,0.00284879,0.00284594,0.00284346,0.00284056,0.00283777,0.00283502,0.0028323,0.00282939,0.00282645,0.00282353,0.00282053,0.00281778,0.00281457,0.00281189,0.00280891,0.00280615,0.00280355,0.00280057,0.00279736,0.00279468,0.00279185,0.00278916,0.00278699,0.0027844,0.00278164,0.00277904,0.00277575,0.00277311,0.00277044,0.00276776,0.00276461,0.00276186,0.00275906,0.00275677,0.00275435,0.00275168,0.00274942,0.00274617,0.00274363,0.00274107,0.00273819,0.0027355,0.00273297,0.00273059,0.00272816,0.00272563,0.00272344,0.00272085,0.00271833,0.00271609,0.00271369,0.00271133,0.00270894,0.0027064,0.00270412,0.00270171,0.0026994,0.00269703,0.00269418,0.00269178,0.00268936,0.00268728,0.00268506,0.00268312,0.00268036,0.00267787,0.00267564,0.00267347,0.00267118,0.00266877,0.00266653,0.00266441,0.00266271,0.00266054,0.00265825,0.0026559,0.00265383,0.00265161,0.00264968,0.00264732,0.00264539,0.00264295,0.00264098,0.00263869,0.00263636,0.00263434,0.00263215,0.00263023,0.00262822,0.00262628,0.00262422,0.00262236,0.00262036,0.00261831,0.00261628,0.00261451,0.00261257,0.00261059,0.00260871,0.00260687,0.00260501,0.00260267,0.00260094,0.00259868,0.00259649,0.00259478,0.00259308,0.00259143,0.00258955,0.00258721,0.00258551,0.00258347,0.00258174,0.00258005,0.00257832,0.00257656,0.00257497,0.00257335,0.00257165,0.00256984,0.00256803,0.00256624,0.00256436,0.00256254,0.0025607,0.00255925,0.00255792,0.00255616,0.00255439,0.0025528,0.00255149,0.00254971,0.00254832,0.00254703,0.00254555,0.00254428,0.00254283,0.00254166,0.00254023,0.00253862,0.00253728,0.002536,0.00253436,0.00253252,0.00253096,0.00252943,0.00252819,0.00252656,0.00252508,0.00252379,0.00252243,0.00252084,0.00251975,0.00251831,0.00251666,0.00251521,0.00251391,0.00251267,0.00251144,0.00251024,0.00250897,0.00250758,0.00250646,0.00250525,0.00250396,0.00250248,0.00250084,0.00249962,0.00249836,0.00249712,0.00249567,0.00249435,0.00249341,0.00249217,0.00249093,0.00248988,0.00248849,0.00248716,0.00248589,0.00248488,0.00248398,0.00248264,0.00248155,0.00248047,0.00247944,0.00247824,0.00247724,0.00247609,0.00247492,0.00247375,0.00247262,0.00247147,0.0024706,0.00246888,0.0024678,0.00246675,0.00246572,0.00246487,0.00246414,0.00246316,0.0024619,0.00246083,0.00245995,0.00245894,0.002458,0.0024569,0.00245599,0.00245492,0.0024538,0.00245275,0.0024515,0.00245035,0.0024494,0.00244835,0.00244731,0.00244625,0.00244542,0.00244473,0.0024439,0.00244295,0.00244172,0.00244076,0.00243982,0.00243907,0.00243812,0.00243725,0.00243639,0.00243555,0.00243464,0.0024338,0.00243277,0.00243145,0.00243048,0.00242965,0.0024285,0.00242736,0.00242654,0.00242575,0.00242483,0.00242391,0.00242298,0.00242193,0.00242078,0.00241981,0.00241862,0.00241771,0.00241663,0.0024155,0.0024146,0.00241367,0.0024125,0.00241116,0.00241008,0.00240902,0.00240729,0.00240548,0.00240419,0.00240253,0.00240081,0.00239908,0.00239741,0.00239517,0.0023929,0.00239052,0.00238798,0.00238505,0.00238189,0.00237904,0.00237565,0.00237201,0.00236829,0.00236374,0.00235923,0.00235403,0.00234858,0.00234288,0.00233677,0.00233038,0.00232438,0.0023171,0.00230981,0.00230198,0.00229329,0.00228476,0.00227551,0.00226594,0.00225598,0.00224564,0.00223392,0.00222138,0.00220881,0.00219607,0.00218362,0.00217038,0.00215639,0.00214182,0.00212751,0.00211102,0.00209334,0.00207535,0.00205704,0.00203711,0.00201595,0.00199504,0.00197257,0.00194934,0.0019241,0.00189856,0.0018717,0.00184359,0.00181463,0.00178446,0.00175283,0.0017206,0.0016859,0.00165096,0.00161563,0.00157919,0.00154094,0.00150137,0.00146059,0.00141945,0.00137825,0.0013367,0.00129363,0.00124908,0.00120422,0.00115974,0.00111459,0.00106967,0.00102485,0.000979928,0.000935386,0.000891343,0.000847616,0.000805325,0.000763099,0.000721605,0.000680876,0.000641258,0.000602415,0.000564516,0.000529212,0.000494813,0.000461726,0.000429703,0.000398541,0.000369522,0.000342268,0.000316742,0.000292816,0.000270237,0.000249087,0.000228998,0.000210083,0.000192775,0.000176194,0.000160682,0.000146541,0.000133475,0.000121171,0.000110027,9.98983e-05,9.04025e-05,8.1898e-05,7.37485e-05,6.64247e-05,5.991e-05,5.34662e-05,4.79254e-05,4.28124e-05,3.80863e-05,3.37871e-05,2.98488e-05,2.64158e-05,2.31778e-05,2.02341e-05,1.77006e-05,1.53314e-05,1.32501e-05,1.13287e-05,9.72139e-06,8.17381e-06,6.84904e-06,5.73097e-06,4.74769e-06,3.91083e-06,3.20076e-06,2.6024e-06,2.08961e-06,1.69223e-06,1.36759e-06,1.08247e-06,8.47588e-07,6.69622e-07,5.17384e-07,4.00954e-07,3.11072e-07,2.34849e-07,1.76165e-07,1.31673e-07,9.73364e-08,7.19841e-08,5.32892e-08,3.90585e-08,2.83402e-08,2.02418e-08,1.4579e-08,1.06446e-08,7.59531e-09,5.37474e-09,3.79538e-09,2.64724e-09,1.84573e-09,1.29566e-09,9.10203e-10,6.14793e-10,4.26199e-10,2.87287e-10,1.936e-10,1.30732e-10,9.00388e-11,6.2542e-11,4.32799e-11,2.90205e-11,1.95416e-11,1.3365e-11,8.69465e-12,5.75237e-12,3.71099e-12,2.54394e-12,1.73591e-12,1.14175e-12,7.46542e-13,4.68409e-13,3.08728e-13,1.94589e-13,1.17053e-13,7.58227e-14,4.59278e-14,2.89681e-14,1.97988e-14,1.19182e-14,7.29947e-15,4.37438e-15  GT:SU:PE:SR ./.:6:6:0

```

And here is the same event as a BEDPE:

```
chr1    1584078 1584538 chr1    1651058 1651572 11  .   -   +   DUP .   SVTYPE=DUP;SVLEN=66541;END=1651069;STRANDS=-+:6;IMPRECISE;CIPOS=-449,10;CIEND=-10...

```

https://github.com/cancerit/BRASS/wiki/BEDPE

The following is a table describing how SVs are represented in BEDPE format:

```
Heading Type    Description
chr1    String  Chromosome of lower coordinate
start1  0-based int Start coordinate of lower coordinate
end1    1-based int End coordinate of lower coordinate
chr2    String  Chromosome of high coordinate
start2  0-based int Start coordinate of high coordinate
end2    1-based int End coordinate of high coordinate
id/name String  ID of event, correlates with VCF
brass_score int Number of aberrant pairs contributing to the rearrangement group.
strand1 [+-]    Strand of end in 'genomic' context - see table
strand2 [+-]    Strand of end in 'genomic' context

```

# Answer
SURVIVOR has a subcommand to convert an SV vcf to bed, although I'm not sure if it's actual bedpe format.

To install:

```
git clone https://github.com/fritzsedlazeck/SURVIVOR.git
cd SURVIVOR/Debug
make

```

Then:

```
./SURVIVOR vcftobed input.vcf $minSize $maxSize output.bed

```

If you run this on your example file, you get:

```
$ ./SURVIVOR vcftobed foo.vcf 1 10000000000000000 bar.bed; cat bar.bed
chr1    1584528 1584528 chr1    1651069 1651069 11  ,   -   +   DUP

```
> 3 votes

---
Tags: vcf, structural-variation
---


# Question
Title: Error: package or namespace load failed for 'ChIPseeker'
How do I solve the problem of loading the bioconductor **Chipseeker** package? I have installed in the following way and I having an error while loading it and I do not know why.

```
BiocManager::install("ChIPseeker")
BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene")
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")

library(ChIPseeker) # produces and error

Error: package or namespace load failed for 'ChIPseeker':
  'TxDb.Hsapiens.UCSC.hg19.knownGene' is not exported from 'namespace: TxDb.Hsapiens.UCSC.hg19.knownGene'

```

Thank you in advance for any help on how to go about this.

```
session info:

R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252
[4] LC_NUMERIC=C                    LC_TIME=German_Germany.1252

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] BiocInstaller_1.32.1                    TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
 [3] TxDb.Hsapiens.UCSC.hg38.knownGene_3.4.0 GenomicFeatures_1.34.8
 [5] AnnotationDbi_1.44.0                    GenomicRanges_1.34.0
 [7] GenomeInfoDb_1.18.2                     IRanges_2.16.0
 [9] S4Vectors_0.20.1                        cellrangerRkit_1.1.0
[11] Rmisc_1.5                               plyr_1.8.4
[13] lattice_0.20-35                         bit64_0.9-7
[15] bit_1.1-14                              ggplot2_3.1.1
[17] RColorBrewer_1.1-2                      Biobase_2.42.0
[19] BiocGenerics_0.28.0                     Matrix_1.2-14

loaded via a namespace (and not attached):
Error in x[["Version"]] : subscript out of bounds
In addition: Warning message:
In FUN(X[[i]], ...) : DESCRIPTION von Paket 'digest' fehlt oder zerstört

```

# Answer
> 1 votes
Simple, use the latest version of R and Bioconductor, now on version 3.6.1 and Bioconductor 3.9. Those are the versions that are supported. If you can't upgrade R, check that `BiocManager::valid()` is true, you might have ended in some inconsistency between Bioconductor versions.

Looking at your `sessionInfo`, you seem to have some problems with the installed libraries, check and correct those too if possible

---
Tags: r, bioconductor, chip-seq
---


title: What do the read colors in IGV mean?

text:
# Question
Title: What do the read colors in IGV mean?
I was looking at a bam file in the IGV viewer and saw:

What do the different read colors mean? Why is one read a light blue color, another green, another aquamarine (?), another purple and another two orange? What do these different read colors represent? The IGV docs have sections on interpreting color by insert size, pair orientation, bisulphide mode, splice junctions and sashimi plot, but none of them seems to be describing what I'm looking at there.

More generally, where are all these default behaviors documented? Are they somewhere in http://software.broadinstitute.org/software/igv/UserGuide and I'm missing them?

# Answer
> 5 votes
Could it be that these are reads with their mates found on other chromosomes? That's the part of the IGV documentation which seems to make sense in this context:

> for paired end reads that are coded by the chromosome on which their mates can be found

---
Tags: sequence-alignment, visualization, igv
---


# Question
Title: Error while calling bcftools mpileup - Failed to open -: unknown file type
I have sequenced a bacterial genome with a GridIon from ONT. Basically what I want to check is whether or not trimming 50 bps at the beginning of the reads will improve alignment against the reference genome and ultimately the call of a consensus sequence.

For that what I have done so far is:

1. concatenated all the fastqs generated from the Gridion into a single (very large) fastq file: `ref494_cat_all.fastq`
2. aligned this concatenated fastq against the reference genome `ref_genbank.fasta` with `minimap2`: `~/Downloads/minimap2-2.17_x64-linux/minimap2 -ax map-ont ref_genbank/ref_genbank.fasta ref494_cat_all.fastq > ref494_aligned_trim50.sam`
3. Performed the general `samtools` steps to covert sam to bam, index and sort

Then, in order to get the consensus sequence I first try to generate a `.vcf.gz` file, according to this tutorial

```
bcftools mpileup -f ref_genbank/ref_genbank.fasta ref494_aligned.sorted.bam | bcftools call -mv -Oz -o calls.vcf.gz

```

Immediatelly after calling the command above I got:

```
Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid
[mpileup] 1 samples in 1 input files

```

Then, after a minute or two, I got:

```
Failed to open -: unknown file type

```

I'm pretty sure that there's no issue with `.bam` file, as I could visualise the alignment with IGV or tablet.

Does anyone know what the issue can be?

# Answer
The `-` is usually used to mean standard input when reading data. So the message would suggest that the second `bcftools` command fails to read from stdin, so it fails to read the output of the first `bcftools` command:

```
$ echo "foo bar" | /ref/software/bcftools-1.8/bcftools call -mv -Oz -o calls.vcf.gz
Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid
Failed to open -: unknown file type

```

So there is something wrong with either your input file or with the output produced when running the first `bcftools` command on that input file. We can't really help more unless you share the file.
> 1 votes

# Answer
It's possible you're looking at instructions for a recent version of samtools / bcftools, whereas your computer has an older version. Can you try `samtools mpileup -g -f ...` to output the pileup as BCF format instead of `bcftools mpileup -f ...`?

Update: it's good to know that the samtools / bcftools pipelines are at least producing consistent results, which suggests that the error is not associated with the output format on the mpileup side of the pipeline, and also good to see that you've identified that this was a memory issue.

The most frequent cause of memory issues I've had with sorted BAM files has been trying to figure out variants in extremely high coverage data (e.g. over 1000X), which would make sense for a bacterial isolate sequenced on a GridION. I use my own script to carry out a reservoir sampling on mapped BAM format files, which reduces coverage down to a more manageable level. Here's an example command to normalise to 100X coverage where possible:

```
samtools view reads.bam | ./samNormalise.pl -coverage 100 | \
   samtools view -b > reads_normalised.bam

```
> 1 votes

---
Tags: sequence-alignment, sam, vcf, samtools, nanopore
---


# Question
Title: Why is 42 the seed used in RunPCA() in Seurat?
The `RunPCA()` in Seurat 3.0 function has `seed.use=42` by default. How does one choose a seed number? Is this chosen simply because it is the Answer to the Ultimate Question of Life, the Universe, and Everything? Or maybe Lewis Carroll's curious obsession with the number?

# Answer
Given that it is virtually impossible for a human to predict what random numbers will be generated given a certain seed, and what effects they will have for a given application, the choice of the seed can be considered arbitrary.

What matters is that fixing the seed should fix the behaviour of the application. What seed has been chosen does not really matter, as long as it is documented.

The choice of `42` could therefore well come from one of the reasons you propose. Only the authors of the programs could confirm.
> 3 votes

---
Tags: r, seurat
---


# Question
Title: Convert BAM to properly paired FASTQ files
I thought I had figured this one out. But apparently not.

I need to convert a BAM file of paired-end alignments to two FASTQ files of paired reads to realign them, with a twist: I only want reads that fall within a defined region.

Currently I am using the following command (the region is an example; it also fails on different regions; at any rate, this is on hs37d5 but presumably the issue would exist for other references too):

```
samtools view -b -u -@ 48 alignments.bam 21:31,831,502-31,896,094 \
| samtools collate -n 128 -u -O -@ 48 - /tmp/bam-to-fastq- \
| samtools fastq -F 0x900 -@ 48 \
    -0 /dev/null -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz -

```

(Apart from the region filtering, this command is essentially the same that word of god tells us to use.)

The resulting files are not properly paired:

```
⟩⟩⟩ diff <(gunzip -c reads_R1.fastq.gz | sed -n -e 1~4p -e 10q) \
1 ⟩   <(gunzip -c reads_R2.fastq.gz | sed -n -e 1~4p -e 10q)
1,2d0
< @A00346:17:H52H7DSXX:3:2142:32316:24518
< @A00346:17:H52H7DSXX:3:1153:6705:26412
3a2,3
> @A00346:17:H52H7DSXX:3:1173:9118:3129
> @A00346:17:H52H7DSXX:3:1109:29414:11757

```

As far as I can see from performing a spot check, these are all reads which are marked as paired in the BAM file, but for which only one mate falls within the filtered region.

And running `bwa mem` on the files gives me the error:

> \[mem\_sam\_pe\] paired reads have different names: \[…\]

How can I ensure that only properly paired reads are included in the FASTQ files?

# Answer
As noted in the comments, the problem is “some reads fall in the target region but their pairs fall outside it”, leading to non-trivial numbers of singleton reads coming out of `samtools collate`.

```
… | samtools fastq -F 0x900 -@ 48 \
    -0 /dev/null -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz -

```

Your `samtools fastq` command is not doing anything to siphon off these singleton reads. You should add `-s /dev/null` to get rid of them.

See the discussion in samtools/samtools#874 (especially the part starting here). I have not looked at how the recent changes from jkbonfield, which will land in an upcoming samtools version (after 1.9), might affect this.
> 6 votes

---
Tags: sam, samtools, fastq
---


title: Raw Data frame manipulation in python

text:
# Question
Title: Raw Data frame manipulation in python
Using python 3 I need to process qPCR sequencing raw data outputs by searching for the first occurrence of a user defined string and then making a new data frame using all lines after that string. I am trying to find solutions in the pandas doc but so far unsuccessful.

This is a raw output .csv file that I need to process. (couldn't paste complete csv as exceeds character limit, this is lines 40-50 and am hoping this text is useful?). I need to tell pandas to create a new data frame that 1. starts at the line containg the first occurance of str("Sample Name") with that line as header and containing all lines following. And then 2., only including columns ("Sample Name"), ("Target Name"), ("CT").

Could someone please help me so that I can use python to help me analyze biological data? Many thanks, Luke

```
40,Quantification Cycle Method,Ct,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

41,Signal Smoothing On,true,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

42,Stage where Melt Analysis is performed,Stage3,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

43,Stage/ Cycle where Ct Analysis is performed,"Stage2, Step2",,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

44,User Name,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

45,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

46,Well,Well Position,Omit,Sample Name,Target Name,Task,Reporter,Quencher,Quantity,Quantity Mean,SE,RQ,RQ Min,RQ Max,CT,Ct Mean,Ct SD,Delta Ct,Delta Ct Mean,Delta Ct SD,Delta Ct SE,Delta Delta Ct,Automatic Ct Threshold,Ct Threshold,Automatic Baseline,Baseline Start,Baseline End,Amp Status,Comments,Cq Conf,CQCONF,HIGHSD,OUTLIERRG,Tm1,Tm2,Tm3,Tm4

47,1,A1,False,WT1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,23.357698440551758,23.4766845703125,0.5336655378341675,,,,,,True,20959.612776965325,True,3,17,Amp,,0.9588544573203085,N,Y,N,81.40960693359375,,,

48,2,A2,False,WT1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,24.05980110168457,23.4766845703125,0.5336655378341675,,,,,,True,20959.612776965325,True,3,15,Amp,,0.9592687354496955,N,Y,N,81.40960693359375,,,

49,3,A3,False,WT1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,23.012556076049805,23.4766845703125,0.5336655378341675,,,,,,True,20959.612776965325,True,3,16,Amp,,0.9592714462250367,N,Y,N,81.40960693359375,,,

50,4,A4,False,fla11fla12-1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,23.803699493408203,24.419523239135742,0.5669151544570923,,,,,,True,20959.612776965325,True,3,17,Amp,,0.9671570584141241,N,Y,N,81.40960693359375,,,

```

This is the code that I have so far:

```
import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

data = pd.read_excel ("2019-02-27_161601 AtWAKL8 different version expressions.xls", sheet_name='Results').fillna(0)

data.to_csv('df1' + '.csv', index=True)

df1 = pd.read_csv ("df1.csv")

```

# Answer
> 2 votes
```
data = pd.read_excel ("2019-02-27_161601 AtWAKL8 different version expressions.xls", sheet_name='Results').fillna(0)

data.to_csv('df1' + '.csv', index=True)

df1 = pd.read_csv ("df1.csv")

```

This doesn't make much sense ... you read in Excel to a dataframe, you then read out to csv - the dataframe doesn't have an index, but you include an index and then read it in again. This could do more harm than good.

Anyway your supervisor is right it is very easy to parse this data using pandas, although there will be more expertise doing this via R.

I don't want to peer through the gloom of a CSV output, but I would like to highlight how easy it is to parse this stuff into a dataframe. The dataframe '`data`' is fine.

1. **Parse** You can simply remove the top rows of bumf `data = data.iloc[6:]` (I might have miss counted, plus see below),
2. **Headers** You will probably have to remove the "headers" row to a list and then introduce back onto the dataframe using `data.columns =
['header1', 'header2' ...]` ... rather this is the list (array in Perl) which you can automatically generate using a magic `to_list` type pandas command which I forget, its not `list` \- its cuter
3. **Columns** You next might use `drop` to remove the columns you don't want (that might be alot of columns) or else make a subset of the dataframe into a new dataframe, you can also use `iloc` for columns too.
4. **Index** You then assign A1, A2, ... etc as an index `data.index = ['A1',
'A2', ]`...etc.. you can reshape column 3 to automate this
5. **to\_numeric** You might have to assign the CT values as floats or better `to_numeric`
6. **Extrapolation** You would then define a subroutine (or whatever they call them in Python) so perform the extrapolation and use the `apply` command on your dataframe and assign it to a new column, where the input is your CT values.
7. **Output** Finally, issue your `to_csv` command `to_pickle` is also a good idea because you have a df that is worth importing at a future date.
8. **Plotting** this is easy using `data.plot` which appears to be some sort of inbuilt matplotlib, keep in mind seaborn DONT try using R's `ggplot2` within Python using rpy2 .... its a NIGHTMARE to get working and doesn't support the latest `ggplot2` (!!!). R s `ggplot2` is stronger than seaborn or matplotlib (providing you are working within R), but these two are good for quick graphs.

If I had time I'd write the code for you, but I'm afraid I'll have to leave it there. I'm still in debt with a Perl efetch script I promised a previous OP.

In summary, Python pandas is fabulous stuff for this purpose, I love it. Doing this via Perl and "subprocessing" inside R... its doable, well in this particular case I suspect its just as easy, but generally its alot more code and often fiddly.

---
Tags: python, data-retrieval, jupyter
---


# Question
Title: R error: could not find function "pathwayAnalysis"
```
## ------------------------------------------------------------------------
library(limma)
library(statmod)
library(GEOquery)
#library(GEOmetadb)
require(Biobase)
library(ggplot2)
library(rio)
library(affy)
library(affyPLM)
library(reshape)
library(stringr)
library(cowplot)
#library(makecdfenv)
#library(aroma.affymetrix)
#library(RTN)
library(gProfileR)
library(corpcor)
library(viper)
library(rio)
library(dplyr)
library(igraph)
library(RColorBrewer)
library(reshape2)

## ------------------------------------------------------------------------
gse89988 <- getGEO(filename = "GSE89988_family.soft.gz", GSEMatrix = TRUE)

## ------------------------------------------------------------------------
# define samples
samples_nt <- c('GSM2394292', 'GSM2394303')
samples_pam <- c('GSM2394295', 'GSM2394306')
samples_lps <- c('GSM2394296', 'GSM2394307')

# make table with samples metadata
samples <- c(samples_nt, samples_pam, samples_lps)
metadata <- as.data.frame(samples)
colnames(metadata) <- 'sample_id'

# add group information both with the group acronym as a dummy variable -- which
# will be useful for the DEA; also add color for differentiate groups in plots
metadata$group <- c(rep(0, 2), rep(1, 2), rep(2, 2))
metadata$group_name <- c(rep('NT', 2), rep('PAM', 2), rep('LPS', 2))
metadata$group_color <- c(rep('gray', 2), rep('tomato', 2), rep('dodgerblue', 2))

# get positions in GSE object
match(samples, names(GSMList(gse89988)))

# with the following I confirm that the samples correspond to the expected samples
# and get that expression values are already normalized and log2-transformed
GSMList(gse89988)[[49]]
GSMList(gse89988)[[60]]
GSMList(gse89988)[[52]]
GSMList(gse89988)[[63]]
GSMList(gse89988)[[53]]
GSMList(gse89988)[[64]]

## ------------------------------------------------------------------------
#First, we need to make sure that all of the GSMs are from the same platform
gsmplatforms <- lapply(GSMList(gse89988),function(x) {Meta(x)$platform_id})
head(gsmplatforms)
# from the output of the command above we see that indeed all GSMs derive from
# the same platform "GPL17543"

# We can filter the original GSMList to include only those GSMs with the GPL96 platform and use this list for further processing
gsmlist = Filter(function(gsm) {Meta(gsm)$platform_id=='GPL17543'},GSMList(gse89988))
length(gsmlist)

# get the probeset ordering
probesets <- Table(GPLList(gse89988)[[1]])$ID

# make the data matrix from the VALUE columns from each GSM
# being careful to match the order of the probesets in the platform
# with those in the GSMs
data.matrix <- do.call('cbind', lapply(gsmlist,function(x)
                                      {tab <- Table(x)
                                      mymatch <- match(probesets,tab$ID_REF)
                                      return(tab$VALUE[mymatch])
                                     }))
data.matrix <- apply(data.matrix, 2, function(x) {as.numeric(as.character(x))})

# go through the necessary steps to make a compliant ExpressionSet
# and filter for the samples to be analyzed
rownames(data.matrix) <- probesets
colnames(data.matrix) <- names(gsmlist)
data.matrix <- data.matrix[, samples]
pdata <- data.frame(samples = samples)
rownames(pdata) <- samples
pheno <- as(pdata,"AnnotatedDataFrame")
eset <- new('ExpressionSet', exprs = data.matrix, phenoData = pheno)

## ------------------------------------------------------------------------
load_gene_annotation <- function() {

  ifile <- 'GPL17543-38936.txt'
  df <- read.delim(ifile, comment.char = '#')

  return(df)

}

gene_annotation <- load_gene_annotation()

## ------------------------------------------------------------------------
apply_dea <- function(eset, fc_cutoff = 2, adj_p_cutoff = 0.05) {

  # output directory
  ODIR <- 'results/step1/'
  dir.create(ODIR, showWarnings = FALSE, recursive = TRUE)

  # make design matrix
  design <- model.matrix(~ 0 + factor(c(0, 0, 1, 1, 2, 2)))
  colnames(design) <- c("group0", "group1", "group2")

  # plot with the normalized array weights
  aWeights <- arrayWeights(eset, design = design)
  ofile <- paste0(ODIR, 'array_weights_gse89988_samples_used_in_dea.pdf')
  pdf(ofile, width = 10, height = 5)
  dotchart(aWeights, labels = samples,
           groups = as.factor(metadata$group_name),
           pch = 19,
           color = metadata$group_color,
           main = "Array Weights")
  dev.off()

  # linear model fit
  fit <- lmFit(eset, design)
  contrast.matrix <- makeContrasts(group1-group0, group2-group0,
                                   levels = design)
  fit2 <- contrasts.fit(fit, contrast.matrix)
  fit2 <- eBayes(fit2)

  # save
  ODIR='deliverables/'
  dir.create(ODIR, showWarnings = FALSE)
  ofile=paste0(ODIR, 'fit_for_figures_7_and_8.rds')
  saveRDS(fit2, ofile)

  # extract results, annotate and visualize DEA results

  extract_annotate_and_viz_dea_results <- function(coef_index, ref, treat) {

    # Extract DEA results for a given comparison
    df <- topTable(fit2, coef = coef_index, adjust = "BH", number = Inf)
    df$reference <- ref
    df$treatment <- treat

    # remove genes with missing values
    df <- df[complete.cases(df), ]

    # annotate with gene information -- gene annotation has many fields which
    # are not relevant here so I will make a reduced gene information version
    fields_included <- c('ID', 'Species', 'Source', 'Entrez_Gene_ID', 'Symbol')
    df <- merge(df, gene_annotation[fields_included], by.x = 0, by.y = 'ID',
                sort = FALSE)
    colnames(df)[1] <- 'probe_id'

    # prepare data for visualization through volcano plot
    cond1 <- df$adj.P.Val < adj_p_cutoff
    cond2d <- df$logFC < -log2(fc_cutoff)
    cond2u <- df$logFC > log2(fc_cutoff)
    log10p <- -log10(df$adj.P.Val)
    down <- nrow(df[cond1 & cond2d, ])
    up <- nrow(df[cond1 & cond2u, ])
    df$category <- ifelse(cond1 & cond2d, 'down-regulated',
                          ifelse(cond1 & cond2u, 'up-regulated', 'non-regulated'))
    abs_max_logfc <- max(abs(df$logFC))

    # volcano plot
    g <- ggplot(df, aes(x = logFC, y = -log10(adj.P.Val), colour = category)) +
          geom_point(alpha = 0.4, size = 1.75) +
          scale_color_manual(values = c('dodgerblue', 'gray', 'tomato')) +
          xlim(-abs_max_logfc, abs_max_logfc) +
          ylim(c(0, max(-log10(df$adj.P.Val)))) +
          xlab("log2 fold change") +
          ylab("-log10(adjusted p-value)") +
          annotate("text", x = -4, y = 3, size =3,
                   label = paste0(down, "\ndown-regulated DEGs"),
                   color = 'dodgerblue') +
          annotate("text", x = 4, y = 3, size = 3,
                   label = paste0(up, "\nup-regulated DEGs"),
                   color = 'tomato') +
          annotate("text", x = -abs_max_logfc, y = 0, size = 2,
                   label = paste0("DEGs defined with:\nAdj. p-value < ",
                                  adj_p_cutoff,
                                  '\nFold change > ',
                                  fc_cutoff),
                   color = 'black', hjust = 0) +
          ggtitle(paste0(ref, ' vs. ', treat)) +
          theme(plot.title = element_text(hjust = 0.5))

    # save plot to file
    ofile <- paste0(ODIR, 'volcano_plot_', tolower(ref), '_vs_', tolower(treat),
                    '.pdf')
    ggsave(ofile, width = 7, height = 5)

    return(df)

  }

  res <- rbind(extract_annotate_and_viz_dea_results(1, 'NT', 'PAM'),
               extract_annotate_and_viz_dea_results(2, 'NT', 'LPS'))

  # export table with results
  ofile <- paste0(ODIR, 'differential_expression_analysis_results_combined.xlsx')
  export(res, ofile)

}

## ------------------------------------------------------------------------
apply_dea(eset)

```

After all these i have to run this

`gsca <- pathwayAnalysis(fit2, pvaluecutoff=0.05)`,where fit2 is made after running `apply_dea(eset)`

As I run the above code i get this

> Error in pathwayAnalysis(fit2, pvaluecutoff = 0.05) : could not find function "pathwayAnalysis"

I would share the .rdata file for reproduciblty
rdata

The gse & the annotation files

I would be really glad if i get to know what am i doing wrong.

# Answer
It seems like you are missing function "pathwayAnalysis" presumably from package BLMA, i.e. `library('BLMA')` should resolve your error once you install that package.

-- edit 1 --

Ha, stupid me. Package `BLMA` contains file "pathwayAnalysis.R", not the function. A package that contains a FUNCTION called `pathwayAnalysis` is IntClust. You should try that one.

-- edit 2 --

Alright, stupid me again, now I have not checked that IntClust has capitalized `PathwayAnalysis`, you need `pathwayAnalysis`. I found one more function matching the name form package `metaX`. However, it could be that it's function defined somewhere in their code. You can try to search the source code for example using grep or you try to use google to find packages that contain the function of this name (this is exactly what I was doing, many you should be a bit more careful than I was).
> 1 votes

---
Tags: r, microarray
---


# Question
Title: Subset a protein fasta based on sequence length?
I have a strange protein fasta file. There are several entries for the same gene and I need to extract the entry with the longest width.

```
>lcl|NW_017095468.1_prot_XP_017786561.1_8 [gene=LOC108569215] [db_xref=GeneID:108569215] [protein=NADH dehydrogenase [ubiquinone] 1 a
lpha subcomplex subunit 9, mitochondrial] [protein_id=XP_017786561.1] [location=complement(join(88030..88986,89044..89243,89412..8944
8))] [gbkey=CDS]
MAAVIFTGAQLLKQQSGLVGIAYIRVNNYSSDAKVYNLASLKRGTGGRSSFNGIVATVFGSTGFIGRYVANKLGKIGTQL
ILPYRGNNYETMRLKLCGDLGQVLFQPYYLKDLESIDKALKYSNVVINLVGRDWETKNFSFQDVHVKGARDLARAAKKAG
VEKFIHLSALNCDGPNEAIFSRQGSKFLSSKWEGEQAVLEEFPEATIFRPSDVYGQEDRFLRYYAHNWRRQGQYMPMWKN
GEATIKQPVHVSDLAAGIVAAIKDPEAAGKVYQAVGPRRYQLNELVDWFYRVMRKDSEWGYKRYDMKYDMFFKLKVSLTQ
KFSPAFPVGNLHWERLEREFVTDRVNQAIPTLEDLGVNLRRMEDQVPWELKPFTYGLYFGTDSEEPVVEPKPPKYVS
>lcl|NW_017095468.1_prot_XP_017770525.1_9 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X1] [protein_id=XP_017770525.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,108885..109058,109127..109237,109573..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD
DGPLSKILCEEVQENEDTQPNMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMSLVSMIDPSPDWIVGVSGLEL
CLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHPIRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEK
NCVAQVDVSEEDGGVEGDKCEMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERECGPSDRRPLQESKECEGRAW
SNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDTEDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDN
SAVVCGFPRDEGGCMSNVDRWYFDVIKGNCDIFSYSGCQGNKNNFKTLERCENVCDSYKKELLANRTAYKRQLGVTVSGV
LSYNLHHMQNDDADNCVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGFSTKHRFIRVHPSNGGKRCPQ
KILKKRKCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCVMTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHR
TEERPCSLLACPE
>lcl|NW_017095468.1_prot_XP_017771299.1_10 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X1] [protein_id=XP_017771299.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,108885..109058,109127..109237,109573..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD
DGPLSKILCEEVQENEDTQPNMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMSLVSMIDPSPDWIVGVSGLEL
CLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHPIRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEK
NCVAQVDVSEEDGGVEGDKCEMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERECGPSDRRPLQESKECEGRAW
SNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDTEDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDN
SAVVCGFPRDEGGCMSNVDRWYFDVIKGNCDIFSYSGCQGNKNNFKTLERCENVCDSYKKELLANRTAYKRQLGVTVSGV
LSYNLHHMQNDDADNCVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGFSTKHRFIRVHPSNGGKRCPQ
KILKKRKCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCVMTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHR
TEERPCSLLACPE
>lcl|NW_017095468.1_prot_XP_017772069.1_11 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X2] [protein_id=XP_017772069.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,109546..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD
DGPLSKILCEEVQENEDTQPNMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMSLVSMIDPSPDWIVGVSGLEL
CLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHPIRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEK
NCVAQVDVSEEDGGVEGDKCEMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERECGPSDRRPLQESKECEGRAW
SNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDTEDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDN
SVTENYFGRIVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGFSTKHRFIRVHPSNGGKRCPQKILKKR
KCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCVMTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHRTEERPC
SLLACPE
>lcl|NW_017095468.1_prot_XP_017772815.1_12 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X3] [protein_id=XP_017772815.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,109573..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD

```

I've tried using BioStrings but I'm unfamiliar with how to index this type of file. Any help appreciated.

# Answer
> 1 votes
Using the `FastaToTbl` and `TblToFasta` scripts I have posted before, you can do:

```
FastaToTbl file.fa |
    awk '{
            len=length($NF);
            if(len > l[$2]){
                l[$2]=len;
                seq[$2]=$0;
            }
          }
          END{
            for(s in seq){
                print seq[s]
            }
        }' | TblToFasta

```

Running this on your example gives:

```
$ FastaToTbl file.fa | awk '{len=length($NF); if(len > l[$2]){l[$2]=len; seq[$2]=$0; }} END{for(s in seq){print seq[s]}}' | TblToFasta
>lcl|NW_017095468.1_prot_XP_017770525.1_9 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X1] [protein_id=XP_017770525.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,108885..109058,109127..109237,109573..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNI
FLRSNGEYQAKNKFKDFLLLVEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVM
QTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTDDGPLSKILCEEVQENEDTQP
NMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMS
LVSMIDPSPDWIVGVSGLELCLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHP
IRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEKNCVAQVDVSEEDGGVEGDKC
EMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERE
CGPSDRRPLQESKECEGRAWSNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDT
EDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDNSAVVCGFPRDEGGCMSNVDR
WYFDVIKGNCDIFSYSGCQGNKNNFKTLERCENVCDSYKKELLANRTAYKRQLGVTVSGV
LSYNLHHMQNDDADNCVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGF
STKHRFIRVHPSNGGKRCPQKILKKRKCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCV
MTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHRTEERPCSLLACPE
>lcl|NW_017095468.1_prot_XP_017786561.1_8 [gene=LOC108569215] [db_xref=GeneID:108569215] [protein=NADH dehydrogenase [ubiquinone] 1 alpha subcomplex subunit 9, mitochondrial] [protein_id=XP_017786561.1] [location=complement(join(88030..88986,89044..89243,89412..89448))] [gbkey=CDS]
MAAVIFTGAQLLKQQSGLVGIAYIRVNNYSSDAKVYNLASLKRGTGGRSSFNGIVATVFG
STGFIGRYVANKLGKIGTQLILPYRGNNYETMRLKLCGDLGQVLFQPYYLKDLESIDKAL
KYSNVVINLVGRDWETKNFSFQDVHVKGARDLARAAKKAGVEKFIHLSALNCDGPNEAIF
SRQGSKFLSSKWEGEQAVLEEFPEATIFRPSDVYGQEDRFLRYYAHNWRRQGQYMPMWKN
GEATIKQPVHVSDLAAGIVAAIKDPEAAGKVYQAVGPRRYQLNELVDWFYRVMRKDSEWG
YKRYDMKYDMFFKLKVSLTQKFSPAFPVGNLHWERLEREFVTDRVNQAIPTLEDLGVNLR
RMEDQVPWELKPFTYGLYFGTDSEEPVVEPKPPKYVS

```

---
Tags: fasta
---


# Question
Title: How to output only the overlapping region of two files?
I have two files of genomic coordinates. I am trying to output only the overlapping region of fileA if it intersects with fileB. Not the original fileA coordinates.

fileA

```
chr1 25 50
chr1 75 200

```

fileB

```
chr1 10 60
chr1 80 90

```

desired output

```
chr1 25 50
chr1 80 90

```

# Answer
bedtools is my go-to program for operations on genomic intervals. In particular, the bedtools intersect operation is what you're looking for here.

```
$ cat fileA
chr1    25      50
chr1    75      200
$ cat fileB
chr1    10      60
chr1    80      90
$ bedtools intersect -a fileA -b fileB
chr1    25      50
chr1    80      90

```

Make sure that the values are separated by tabs rather than spaces, or bedtools will have trouble reading the files.
> 3 votes

---
Tags: genome, interval
---


title: how to remove range from fasta header

text:
# Question
Title: how to remove range from fasta header
Could you please suggest me how I can remove range from fasta header:like these number from below sequences which has some colons indicating the range of the genome

:147010-147657

:149201-149845

```
> FVJT01000004.1:147010-147657 Mycobacterium abscessus subsp. massiliense strain 974
GTGACCAGCAGCGCAGACCGATCCCCCGTCGAGGCGATGGTGAGCCACGCCGAGGGCGCCATCTCCGAGG
ATGACATCGTCGCCGCCGCCCGGGAACGTGCCGTTGACCTGGGCGCGGGCCCGGTGACACCCGCGGTAGG

>FSAX01000002.1:149201-149845 Mycobacterium abscessus subsp. bolletii strain 286 genome
ATGACGACGCCCGTCAGCTTCGATAACAGATTCTTCTCCGCAACGTGGTCTGCTATCGCGCGGCGCGAAC
CTCAGGAGATCAAGGAGCTGCGTGCGGAGAATCTACGCGGCCTGACCGGCCACGTTCTGGAGGTGGGAGC

```

# Answer
> 6 votes
Invest some time in learning unix tools such as `grep`, `sed` and `awk`.

```
sed -r '/^>/s/:[0-9]+-[0-9]+//' input.fasta > output.fasta

```

This line of code matches the pattern `:<numbers>-<numbers>` on lines that start with `>` (FASTA sequence headers) and replaces that pattern with `''`, essentially removing text that matches the pattern.

# Answer
> 1 votes
Alternatively, Perl provides a slightly simpler regular expression syntax (`\d` instead of `[0-9]`) that I prefer over `sed`.

```
perl -pe 's/:\d+-\d+//' input.fasta > output.fasta

```

---
Tags: bash
---


# Question
Title: How to combine multiple files into one file?
I have multiple files (n=86000) with one column each and I want to combine them all into one file with 86000 columns.

I tried the following command

```
Paste *.txt > mean.txt

```

And as a for loop

```
for i in *.txt, do paste $i > mean.txt; done

```

This did not appear to work because there are 86000 to combine, I then divided my files into multiple folders with each folder not containing more than 15000 files but this did not working

Format of my single file is as follows (It has 16000 lines)

```
0.047202943506056
0.0383549703366668
0.0159162908608477
0.0238048058550539
0.0694728928685859
0.0389052658410832

```

What is the bug in the code?

# Answer
> 4 votes
The bug in your loop appears be "`> filename`" per column, this over-rides the last iteration resulting in a file with a single column, which has been over-written 85999 times. Thus the file output you observe should be the single column present in the last file of your array. If you replaced "`> filename`" with "`>> filename`" within the loop you would then append to the file a new column with each iteration. If the loop fell over at some point you could then simply use one of the outputs of `wc filename` to assess how many columns you had (you might have to transponse the data for wc to work (lots of scripts for this in Perl), although there will almost certainly be an easy solution for counting columns in awk.

Personally, I would use a dataframe for this operation (possibly R, definately Python pandas), because 86000 columns is a lot of columns in any other walk of life, even if it is just post-processing. Generally, if you just want bash, awk is a better solution in this case (I don't code with it, but some here certainly do). However, I do get the idea that `paste` will work if each column is an identical length.

BTW `paste` must be lower case, you've got `Paste` in the first line of your code, but its okay in the second bit of code.

BTW BTW, I did not know shell would take \*.txt directly as an array, with Perl/Python you need to "glob" it and assign it as an array/list.

---
Tags: sequence-alignment, phylogenetics, phylogeny, bash
---


# Question
Title: Strange p-value histogram for differential gene expression analysis
I'm trying to perform pretty standard differential expression analysis using RNA-seq data.

I've used Kallisto to perform RNA quantification and am using Sleuth to perform the differential expression analysis.

I have 2 conditions that I'd like to compare, which is stored in the `metadata`:

```
> metadata
    sample cell replicate batch                   path
1:  Base-1 Base         1   1   /path/to/Batch1/Base-1
2:  Base-2 Base         2   2   /path/to/Batch2/Base-2
3:  Base-3 Base         3   2   /path/to/Batch2/Base-3
4:  Test-1 Test         1   2   /path/to/Batch2/Test-1
5:  Test-2 Test         2   2   /path/to/Batch2/Test-2
6:  Test-3 Test         3   3   /path/to/Batch3/Test-3

```

I'm attempting to test for differentially expressed genes between the `Base` and `Test` cells while controlling for `batch` (the sequencing batch). The R code to do this is adapted from Sleuth's manual:

```
sleuth_obj = sleuth_prep(
    metadata,
    extra_bootstrap_summary = TRUE,
    num_cores = 1,
    target_mapping = transcripts_genes, #ENSEMBL information mapping transcripts to genes
    aggregation_column = "ens_gene"
)

# smooth raw Kallisto abundances
sleuth_obj = sleuth_fit(sleuth_obj, ~cell + batch, "full")

# fitting error measurements on reduced model
sleuth_obj = sleuth_fit(sleuth_obj, ~batch, "reduced")

# compare models to calculate differentially expressed transcripts
sleuth_obj = sleuth_lrt(sleuth_obj, "reduced", "full")

# extract results
sleuth_table = as.data.table(sleuth_results(
    sleuth_obj,
    'reduced:full',
    'lrt',
    show_all = FALSE
))
sleuth_table_tx = as.data.table(sleuth_results(
    sleuth_obj,
    'reduced:full',
    'lrt',
    show_all = FALSE,
    pval_aggregate = FALSE
))

```

To check the quality of the null hypothesis, I've checked the p-value histograms of the gene-level and transcript-level analyses, and I get the following:

Neither of these plots show well-behaved p-values, which makes me question the validity of the reduced model. How can I go about adjusting my analysis to get more reliable results?

# Edit following from @swbarnes2's answer:

I've used `fdrtool` to create an empirical null model of Sleuth's likelihood ratio test.

```
library(fdrtool)
local_fdr_txs = fdrtool(so$tests$lrt$`reduced:full`$test_stat, statistic="normal")

```

Which gives the following plots:

Plotting the distribution of these p-values produces:

This still has a hump, but looks more anti-conservative than before.

My remaining question is about whether I can use `fdrtool` to model the null hypothesis as normal, since I'm not certain about the distribution for Sleuth's likelihood ratio test

# Answer
> 1 votes
I found that using the Wald test instead of the Likelihood Ratio Test produced very different results with respect to the p-values.

```
sleuth_obj = sleuth_prep(
    metadata,
    extra_bootstrap_summary = TRUE,
    read_bootstrap_tpm = TRUE,
    num_cores = 1,
    target_mapping = transcripts_genes,
    aggregation_column = "ens_gene"
)

# same model of covariates
sleuth_obj = sleuth_fit(sleuth_obj, ~cell + batch, "full")

# using Wald test instead of LRT
sleuth_obj = sleuth_wt(sleuth_obj, "cellLT1", "full")

# extract results
sleuth_table_genes <- as.data.table(sleuth_results(sleuth_obj, 'cellLT1', 'wt', show_all = FALSE, pval_aggregate = TRUE))
sleuth_table_txs <- as.data.table(sleuth_results(sleuth_obj, 'cellLT1', 'wt', show_all = FALSE, pval_aggregate = FALSE))

```

The p-values are different, but surprisingly the rankings of the genes aren't too different from each other.

# Answer
> 1 votes
Try the protocol here for generating a different (smaller) dispersion estimate

https://www.huber.embl.de/users/klaus/Teaching/DESeq2Predoc2014.html#inspection-and-correction-of-pvalues

---
Tags: rna-seq, differential-expression, kallisto, p-values
---


# Question
Title: CRISPR Sequence Finder and Database Download
I am searching for tools to pull CRISPR Spacers from Bacterial Genomes. I am aware of the CRISPRDB and the corresponding identification tool on the web server.

* Are there other tools for finding CRISPR Spacers?
* Is there a way to download the spacer sequences from CRISPRDB?
* Is there another database that contains all validated / inferred bacterial CRISPR sequences?

# Answer
CRT (CRISPR Recognition Tool) http://www.room220.com/crt/

Crisprs Finder Online Tool https://crispr.i2bc.paris-saclay.fr/Server/

Piler-CR http://www.drive5.com/pilercr/

(Omic Tools List) https://omictools.com/crispr-detection-category
> 2 votes

---
Tags: database, software-recommendation, data-retrieval, crispr, bacteria
---


# Question
Title: How to run 'join' command for multiple files (of 2 types) in a folder
I have two types of files for a gene in a folder

File 1: FOS.tf.txt

```
ABL2
ACTN4
ADGRE5
ADIPOR2
ADRB2
ERCC4
EZR
FAS
FMO4

```

File 2: FOS.tt.txt

```
ABL2
CASP9
CCK
DDIT3
ERCC4
EZR
FAS
FMO4

```

After running command

```
join FOS.tf.txt FOS.tt.txt > FOS_final.txt

```

I got output file: FOS\_final.txt

```
ABL2
ERCC4
EZR
FAS
FMO4

```

My question is that I have these two type of files (tf and tt) for 605 genes in a folder and I want to run join command together for all my genes and get 605 output files.

# Answer
> 1 votes
If all of your `.tf.txt` and `.tt.txt` files are in the same directory, the following command should work.

```
for tffile in $(ls *.tf.txt)
do
    prefix=$(basename $tffile .tf.txt)
    join ${prefix}.tf.txt ${prefix}.tt.txt > ${prefix}_final.txt
done

```

---
Tags: bash
---


title: Protein SCA and DCA

text:
# Question
Title: Protein SCA and DCA
Could you please suggest suitable R packages for statistical coupling and direct coupling analysis of protein sequences?
I know there are Python packages for those: SCA and DCA. I was looking for similar packages in R so that I can continue with my present R workflow.

# Answer
While my suggestion doesn't technically fit your qualifications, I do think it may be useful. I bring this up because you mention you want to continue with your R workflow. I am assuming you don't have anything against Python, but just want to stick with R to keep a tidy project.

The R package `reticulate` (homepage) allows for a very easy interface between Python and R. I use it to run the Python package `scrublet` to remove doublets from scRNA-Seq data before analysis with a popular R package.

Just for reference, here is what using a Python package in R might look like:

```
## setup
Sys.setenv("RETICULATE_PYTHON"="/n/apps/CentOS7/install/pyenv-1.0.0/pyenv/shims/python3")
use_python("/n/apps/CentOS7/install/pyenv-1.0.0/pyenv/shims/python3", required=TRUE)
reticulate::py_config()
scr <- import(module="scrublet", convert=FALSE)

## usage
E150_scrub <- scr$Scrublet(r_to_py(seurat_object@assays$RNA@counts)$T$tocsc(), expected_doublet_rate=0.06)
doublet_results <- E150_scrub$scrub_doublets(min_counts=2, min_cells=3, min_gene_variability_pctl=85)

```

The Python package function calls look like `scr$Scrublet()`. There are some very convenient features like easy conversion of numpy tables to dataframes and the like. If you use RMarkdown, you can even display inline plots from `matplotlib`. I highly recommend giving it a try, it smoothly bridges the gap between the spheres of development of bioinformatic tools.
> 3 votes

---
Tags: r, sequence-analysis, sequence-homology
---


# Question
Title: How can I transform a mapped BAM file into an unmapped BAM file?
In order to use MergeBamAlignment (Picard), I need an unmapped BAM file and a mapped BAM file. I have two mapped BAM files:

* one with reads mapped to the reference I want but without metadata such as read group
* a second one with reads mapped to the incorrect reference but with metadata

When I try to use MergeBamAlignment with these two files, I get errors, such as "Second read from pair not found in unmapped bam". So I am trying to unmap the reads from my second BAM file. How can I do that? Is there a tool for it? Otherwise, I'll have to manually change the bitwise FLAGs, reference names, CIGARS, etc.

# Answer
> 2 votes
Looks like you need Picard's `RevertSam` operation (described in context here, manual here).

```
java -jar picard.jar RevertSam \
     I=input.bam \
     O=reverted.bam

```

# Answer
> 0 votes
What happens if you use the bam with the metadata as the input to your aligner? Maybe it will keep all those fields.

---
Tags: sam, read-mapping, software-recommendation
---


# Question
Title: Pandas automatically rounds GWAS P-value
I am working with a specific GWAS. If I were to run this on the command line

```
grep <rs_id_of_interest> GWAS.txt

```

I would see the GWAS p-value to be on the order of 10^-200. Very tiny p-value.

However, if I were to do the same but in pandas, with the following code,

```
gwas = pd.read_csv(GWAS_file,sep='\t')
gwas_rs_id = gwas.loc[gwas['marker'] == 'rs_id'].copy()

```

and then print gwas\_rs\_id, it shows the gwas p-value to be 0.0. Is there a way to read the file and not have it automatically round the tiny tiny p-values?

Thanks!

# Answer
The solution you likely want is here,

```
pd.set_eng_float_format(accuracy=x, use_eng_prefix=True)

```

x = whatever is required The function set\_eng\_float\_format has been moved around a bit and is now a top level function

---

You might be dealing with maximum likelihood, so convert to log likelihood -200 that sort of thing

Best idea

```
gwas_rs_id['logged'] = np.log(gwas_rs_id.marker)

```

First idea,

```
import decimal as D
....
gwas_rs_id['marker'] = gwas_rs_id['marker'].astype(D.decimal)

```

Last attempt,

```
pd.set_option('display.precision', x) # x whatever dp you need, this might work in combination with a log likelihood

```

Try a `gwas_rs_id.dtype()` check whether you're using float32, float64 or still in object. It should convert to scientific notion (E values) automatically.
> 2 votes

---
Tags: python, phylogenetics, phylogeny, gwas
---


# Question
Title: Using `install.packages` with conda-managed R
I have R installed and managed by conda (miniconda) on my MacBook Pro. The version of R I use most frequently (3.5.1) is installed on the base environment and I have other version-specific environments as well. I did this as it was easier than managing R on homebrew.

The downside I now face is that I cannot use `install.packages()`, `BiocManager::install()` or `biocLite()` to install most packages. `devtools::install_github()` seems to work fine for the few packages I have installed from GitHub.

I have to use `conda install -c [conda-forge|bioconda] [r|bioconductor]-package_name` to install my packages and it breaks my workflow and I need to step outside of R to install packages.

Am I doing anything wrong, or is this the only way packages can be managed when R in installed through conda?

A few sample packages that I needed to install using `conda`:

* tidyverse
* tximport
* deseq2
* apeglm
* maftools
* sva

# Answer
Installing packages outside of an R session is arguably the point of using `conda`. Package management isn't something that should be done within scripts, only the loading of packages. Doing so ensures that your environment setup is somewhat robust, which should be set before you perform actual calculations.

Using `conda install ...` is the best way to ensure environment stability here, and, thankfully, lots of people at Bioconductor have worked very hard to make sure `bioconda` packages fit with the conda framework.

If you need to find packages, you can search online at https://anaconda.org, or from the command line via `conda search <pkg>`.
> 3 votes

---
Tags: r, bioconductor, conda
---


# Question
Title: for loop cutadapt for files on a single directory
I have a bigfolder where i have a lot of fastq.gz files and I want to remove the adapters from all of them.

I am trying then the following loop:

```
for file in LotsofRaws\/cannam/*_R1.fastq.gz; do cutadapt
     -a runs/new/Adapters/IlluNexteraPE.fa
        -o LotsofRaws/RemovedAdapter/Adapted$file
           $file done;  done

```

However I am getting the following error:

```
cutadapt: error:

*It appears you want to trim paired-end data because you provided two input files, but then you also need to provide two output files (with -o and -p) or use the --interleaved option.*

```

But I just want to take the adapters from the files I am not looking to do it in the same way I would do with a single file is there any way to do this because if I try to do `--interleaved` I get:

```
cutadapt: error: Character '/' in adapter sequence 'RTNS/NEW/ADAPTERS/ILLTNEXTERAPE.FA' is not a valid IUPAC code. Use only characters XACGTURYSWKMBDHVN.

```

This is cutadapt 2.3 with Python 3.6.7

and if instead of using the file I directly use the sequences:

```
for file in LotsofRaws/cannam/*;
         do cutadapt -a CTGTCTCTTATACACATCT
               -a AGATCGGAAGAGC
               --interleaved
               -o LotsofRaws/RemovedAdapter/$file   $file ; done

```

I get:

```
FileNotFoundError: [Errno 2] No such file or directory: 'LotsofRaws/RemovedAdapter/LotsofRaws/cannam/NAM7_Glenlea_R1.fastq.gz'
During handling of the above exception, another exception occurred:

```

many times

What should I then do for general variant calling it is appropriate to do the --interleaved option?

Should I just the do the trimming manually or can I loop it in any way that still resembles doing it one by one but without me having to go and write line and line again.

So the files on the folder look like:

```
NAM11_BW415_R1.fastq.gz         NAM13_BW968_R2.fastq.gz     NAM17_AC_Barrie_R1.fastq.gz  NAM3_PI648600_R2.fastq.gz  NAM7_Glenlea_R1.fastq.gz   PI648598_R2.fastq.gz
NAM11_BW415_R2.fastq.gz         NAM13_Frontana_R1.fastq.gz  NAM17_AC_Barrie_R2.fastq.gz  NAM3_PI648842_R1.fastq.gz  NAM7_Glenlea_R2.fastq.gz   PI648611_R1.fastq.gz
NAM11_Reeder_R1.fastq.gz        NAM13_Frontana_R2.fastq.gz  NAM17_BW388_R1.fastq.gz      NAM3_PI648842_R2.fastq.gz  NAM7_Laura_R1.fastq.gz     PI648611_R2.fastq.gz
NAM11_Reeder_R2.fastq.gz        NAM14_BW874_R1.fastq.gz     NAM17_BW388_R2.fastq.gz      NAM4_PI613279_R1.fastq.gz  NAM7_Laura_R2.fastq.gz     PI648672_R1.fastq.gz
NAM11_Webill_R1.fastq.gz        NAM14_BW874_R2.fastq.gz     NAM17_RedFife_R1.fastq.gz    NAM4_PI613279_R2.fastq.gz  NAM7_Marquis_R1.fastq.gz   PI648672_R2.fastq.gz
NAM11_Webill_R2.fastq.gz        NAM14_Glenn_R1.fastq.gz     NAM17_RedFife_R2.fastq.gz    NAM4_PI613280_R1.fastq.gz  NAM7_Marquis_R2.fastq.gz   PI648724_R1.fastq.gz
NAM12_BW431_L002_R1.fastq.gz    NAM14_Glenn_R2.fastq.gz     NAM1_PI648517_R1.fastq.gz    NAM4_PI613280_R2.fastq.gz  NAM8_FL62R1_R1.fastq.gz    PI648724_R2.fastq.gz
NAM12_BW431_L002_R2.fastq.gz    NAM14_Unity_R1.fastq.gz     NAM1_PI648517_R2.fastq.gz    NAM4_PI648524_R1.fastq.gz  NAM8_FL62R1_R2.fastq.gz    PI648768_R1.fastq.gz
NAM12_BW431_L008_R1.fastq.gz    NAM14_Unity_R2.fastq.gz     NAM1_PI648519_R1.fastq.gz    NAM4_PI648524_R2.fastq.gz  NAM8_Muchmore_R1.fastq.gz  PI648768_R2.fastq.gz
NAM12_BW431_L008_R2.fastq.gz    NAM15_BW881_R1.fastq.gz     NAM1_PI648519_R2.fastq.gz    NAM5_PI648523_R1.fastq.gz  NAM8_Muchmore_R2.fastq.gz  PI648860_R1.fastq.gz
NAM12_Penhold_L002_R1.fastq.gz  NAM15_BW881_R2.fastq.gz     NAM1_PI648541_R1.fastq.gz    NAM5_PI648523_R2.fastq.gz  NAM8_Stettler_R1.fastq.gz  PI648860_R2.fastq.gz
NAM12_Penhold_L002_R2.fastq.gz  NAM15_Lillian_R1.fastq.gz   NAM1_PI648541_R2.fastq.gz    NAM5_PI648585_R1.fastq.gz  NAM8_Stettler_R2.fastq.gz  PI648862_R1.fastq.gz
NAM12_Penhold_L008_R1.fastq.gz  NAM15_Lillian_R2.fastq.gz   NAM2_PI648493_R1.fastq.gz    NAM5_PI648585_R2.fastq.gz  PI277012_R1.fastq.gz       PI648862_R2.fastq.gz
NAM12_Penhold_L008_R2.fastq.gz  NAM15_PT245_R1.fastq.gz     NAM2_PI648493_R2.fastq.gz    NAM5_PI648622_R1.fastq.gz  PI277012_R2.fastq.gz       PI648864_R1.fastq.gz
NAM12_Sumai3_L002_R1.fastq.gz   NAM15_PT245_R2.fastq.gz     NAM2_PI648602_R1.fastq.gz    NAM5_PI648622_R2.fastq.gz  PI648478_R1.fastq.gz       PI648864_R2.fastq.gz
NAM12_Sumai3_L002_R2.fastq.gz   NAM16_BA51B92_R1.fastq.gz   NAM2_PI648602_R2.fastq.gz    NAM6_Neepawa_R1.fastq.gz   PI648478_R2.fastq.gz       run_701_150904_7001410_0194_BC79VNANXX_NAM.list
NAM12_Sumai3_L008_R1.fastq.gz   NAM16_BA51B92_R2.fastq.gz   NAM2_PI648609_R1.fastq.gz    NAM6_Neepawa_R2.fastq.gz   PI648487_R1.fastq.gz       run_760_160912_7001343F_0065_AC9EV0ANXX_NAM.list
NAM12_Sumai3_L008_R2.fastq.gz   NAM16_BW928_R1.fastq.gz     NAM2_PI648609_R2.fastq.gz    NAM6_PI648532_R1.fastq.gz  PI648487_R2.fastq.gz
NAM13_BW462_R1.fastq.gz         NAM16_BW928_R2.fastq.gz     NAM3_PI648511_R1.fastq.gz    NAM6_PI648532_R2.fastq.gz  PI648566_R1.fastq.gz
NAM13_BW462_R2.fastq.gz         NAM16_BW963_R1.fastq.gz     NAM3_PI648511_R2.fastq.gz    NAM6_Thatcher_R1.fastq.gz  PI648566_R2.fastq.gz
NAM13_BW968_R1.fastq.gz         NAM16_BW963_R2.fastq.gz     NAM3_PI648600_R1.fastq.gz    NAM6_Thatcher_R2.fastq.gz  PI648598_R1.fastq.gz

```

# Answer
You have various different issues here. First, the first version of the command you show would give a syntax error, not the error you show, so I am guessing what you actually ran was this:

```
for file in LotsofRaws\/cannam/*_R1.fastq.gz; do
    cutadapt -a runs/new/Adapters/IlluNexteraPE.fa \
             -o LotsofRaws/RemovedAdapter/Adapted$file $file done;
done

```

The error there is that you have the word `done` after `$file`. So what `cutadapt` sees is:

```
cutadapt [OPTIONS] $file done

```

So it reads that as being given two input files instead of one. What you wanted to do was this:

```
for file in LotsofRaws\/cannam/*_R1.fastq.gz; do
    cutadapt -a runs/new/Adapters/IlluNexteraPE.fa \
             -o LotsofRaws/RemovedAdapter/Adapted$file $file
done

```

Or, on one line, this:

```
for file in LotsofRaws\/cannam/*_R1.fastq.gz; do cutadapt -a runs/new/Adapters/IlluNexteraPE.fa -o LotsofRaws/RemovedAdapter/Adapted$file $file; done

```

One `done`, not two and separated from the command by a newline or a `;`.

Your next attempt, with `--interleaved` failed because the `-a` option expects an actual sequence, not a file unless you use it with `file:`. See `cutadapt --help` (emphasis mine):

> Finding adapters:
> Parameters -a, -g, -b specify adapters to be removed from each read (or from the first read in a pair if data is paired). If specified multiple times, only the best matching adapter is trimmed (but see the --times option). **When the special notation 'file:FILE' is used, adapter sequences are read from the given FASTA file**.

So, the `-a` expects to be given a sequence unless you use it as `-a file:/path/to/file`. Because you used it as `-a runs/new/Adapters/IlluNexteraPE.fa` instead of `-a file:withruns/new/Adapters/IlluNexteraPE.fa`, it tried to interpret `runs/new/Adapters/IlluNexteraPE.fa` as a sequence and complained because the character `/` isn't valid in a DNA sequence.

Finally, your last lop failed because you're giving the wrong name. The first rule of debugging is "print all the things". When you see a command that doesn't work as you expect, your first step should always be to print out the command and see what is actually being run. When you do a `for file in /long/path/to/file` loop, the entire path is included in your `$file` variable. Try running this:

```
for file in LotsofRaws/cannam/*; do
    echo "cutadapt -a CTGTCTCTTATACACATCT -a AGATCGGAAGAGC --interleaved -o LotsofRaws/RemovedAdapter/$file $file";
done

```

You will see that that will print commands like this:

```
cutadapt -a CTGTCTCTTATACACATCT -a AGATCGGAAGAGC --interleaved \
         -o LotsofRaws/RemovedAdapter/LotsofRaws/cannam/NAM11_BW415_R1.fastq.gz \
          LotsofRaws/cannam/NAM11_BW415_R1.fastq.gz

```

What you need is to take the name of the file, without the path. And you don't want to escape the slash after `LotsofRaws`. So putting all this together and going back to your first attempt, what you want is something like this:

```
for f in LotsofRaws\/cannam/*_R1.fastq.gz; do
    fname=$(basename file)
    cutadapt -a "file:runs/new/Adapters/IlluNexteraPE.fa"
             -o "LotsofRaws/RemovedAdapter/Adapted${fname}" $file
done

```
> 0 votes

---
Tags: shell, cutadapt
---


# Question
Title: How are UMIs used to dedupulicate in Drop-seq tools?
The module `DigitalExpression` which is part of the popular Drop-seq tools digitally count gene transcripts.

The manual is not very clear on how exactly it resolves the UMI information. Maybe I'm missing something but I was wondering if UMI are used in isolation or in conjunction with the barcode and or transcript?

The manual i.e. Drop-seq Alignment Cookbook only states:

```
To digitally count gene transcripts, a list of UMIs in each gene, within each cell, is assembled, and UMIs
within edit distance = 1 are merged together. The total number of unique UMI sequences is counted,
and this number is reported as the number of transcripts of that gene for a given cell.

```

To clarify my point, I'll use an example. Say my reads are only 4 bp and read 1 has the cell barcode at position 1+2 and UMI at position 3+4. read 2 constitutes the transcript read.

Given the following example reads:

```
seqence1:   read1 = AATT  read2 = ATGC
seqence2:   read1 = AATT  read2 = TTTT

```

UMI and barcode are the same but read2 is different, would they be collapsed to 1 UMI count or 2?

Similarly as case with same UMI and transcript but different cell barcode:

```
seqence1:   read1 = AATT  read2 = ATGC
seqence2:   read1 = GGTT  read2 = ATGC

```

Also of note, UMI deduplication is the default setting as of the parameter below:

```
OUTPUT_READS_INSTEAD=Boolean  Output number of reads instead of number of unique molecular barcodes.  Default value:
                              false. This option can be set to 'null' to clear the default value. Possible values:
                              {true, false}

```

# Answer
> 3 votes
> UMI and barcode are the same but read2 is different, would they be collapsed to 1 UMI count or 2?

If the two reads share cell barcodes and UMI sequences and align to the same **gene**, (even if they are different sequences) then they originate from the same RNA molecule, and should only count as 1.

https://support.10xgenomics.com/single-cell-gene-expression/library-prep/doc/user-guide-chromium-single-cell-3-reagent-kits-user-guide-v3-chemistry

# Answer
> 1 votes
The code does EXACTLY what it says in the first quote.

Let's say you have 5 reads on a cell barcode that align to gene A. Those reads have the following molecular barcodes:

AAAA,BBBB,AAAA,CCCC,DDDD.

There are 4 UNIQUE molecular barcodes.

The only complication to this is that we collapse barcodes with an edit distance of 1, which takes care of over-inflation of UMI counts due to PCR/Sequencing errors. If we added a barcode with the sequence AAAB to the original list, it would be collapsed into AAAA as it has an edit distance of 1, and you'd still have 4 unique molecular barcodes.

---
Tags: scrnaseq, umi
---


# Question
Title: What are the right parameters to trim a small RNA transcriptome with trimmomatic?
I'm having some problems with finding the right parameters to trim my small RNA Illumina reads (51 nt long) with Trimmomatic.

Before trimming, one of the samples (21M reads) looks like this:

So for my understanding, it's quite good and the last two issues can be easily solved with adapter trimming.

I then tried to trim it with the command

```
java -jar /opt/Trimmomatic-0.38/trimmomatic-0.38.jar SE -threads 24 -phred33 FemaleMito3.fastq FemaleMito3_trimmed.fastq ILLUMINACLIP:adapters.ultimate.fa:2:30:10 AVGQUAL:25 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:18

```

but I got the following result:

So now I have a problem with the sequence quality per tile (why?), an overall drop of the quality at the end of the reads and a warning with the sequence length distribution.

Then I tried trimming just the adapters with this command:

```
java -jar /opt/Trimmomatic-0.38/trimmomatic-0.38.jar SE -threads 24 -phred33 FemaleMito1.fastq FemaleMito1_noadapters.fastq ILLUMINACLIP:adapters.ultimate.fa:2:30:10 AVGQUAL:25

```

but again I get some new issues:

Again the per tile sequence quality issue, and now a warning about the GC content.

I didn't find anywhere some precise instruction on how to trim a small RNA library with Trimmomatic. **How should I properly trim my small RNA library with Trimmomatic?**

# Answer
The command where you trim with adapters and by quality is perfectly fine. That FastQC isn't perfectly happy is expected. It's a tool made for whole genome sequencing QC, you SHOULD see a number of "fail"s with any RNA-seq protocol. Steps that should always fail in RNA-seq:

* per-base sequence content (there's "random priming" performed that isn't completely random, so the first ~8 bases show a typical bias pattern).
* per-sequence GC content (particularly in sRNA-seq you're enriching for a subset of genes, they won't have some nice GC distribution and that's perfectly OK).
* Sequence duplication level (if this doesn't fail then everything is uniformly expressed...which would be quite unlikely)
* Overrepresented sequences (as above, though this will be worse for sRNA-seq since reads will often cover the entire transcript, so there's even less sequence diversity)
> 3 votes

---
Tags: rna-seq, transcriptome, trimming, data-preprocessing
---


title: How to simulate nanopore reads?

text:
# Question
Title: How to simulate nanopore reads?
I have looked already here: Tools for simulating Oxford Nanopore reads . This doesn't answer my question, because it lists a few Nanopore read simulators, but I have specific problems with each of them.

I am trying to simulate some Nanopore reads. I could use NanoSim, but that tool seems to be very difficult to use. I only have a reference genome, but you need training data as well? anyone have any experience with this? I can only find bacterial and yeast data (and I want to try this on *Arabidopsis thaliana* reference)

Another promising tool I found, was a paper for SNaReSim, reference here: https://ieeexplore.ieee.org/document/8031171. However I cannot seem to find the tool itself, which is also mentioned in the other stackexchange post I mentioned above.

Does anyone have experience simulating PacBio reads with NanoSim, SNaReSim or maybe another tool? Because I am not sure about which one to use. Others like SiLiCO or ReadSim were updated 2 years ago and I believe that Nanopore technology has improved a lot in that time.

To clarify a bit on why I want to simulate data: I just want to test several tools about SV calling and simulating for a school project. All of the points mentioned in the comments I will use for the discussion, but I wanted to focus on 1 read simulator per technology and didn't know what to use, because Nanopore is difficult to simulate in that sense. That was already known to me.

# Answer
It is important to train an error model as NanoSim does, as we do not fully understand the error processes involved in both the nanopore sequencing process and the basecalling process. Any sort of read simulator that does not use a species-specific model is just not going to produce realistic simulated reads.

As far as datasets go, I would recommend looking at this assembly and the associated data uploaded on EMBL:

> Michael, Todd P., et al. "High contiguity Arabidopsis thaliana genome assembly with a single nanopore flow cell." Nature communications 9.1 (2018): 541.

There are several data types there, as well as an assembly to train the models on.

## Is it OK to train a model using reads from one individual and a reference assembly from another individual?

Generally, it is not good to map reads from one individual to the reference genome of another individual since this introduces something called reference bias. When we align reads to a reference genome generated from the same reads, or just from other reads from the same individual we remove this problem. All true reads will map to correctly assembled parts of the genome.

## Why can we train a model for nanopore reads on a reference generated from ... nanopore reads? Won't we get a bad model for indels?

Yes, well, maybe. If the reference genome has been corrected to remove indels using a tool like pilon, then you can generally presume that the reference genome has very few to no indels present as a result of sequencing error from the Nanopore/PacBio reads. As a result, the model derived from aligning reads should have the correct error profile for indels.

Again, it is important that the long-read-assembly was corrected with Illumina shotgun reads from the same individual from which the long reads were derived. Otherwise the corrected reference would have a chimaera of indels and SNPs from multiple individuals.
> 5 votes

# Answer
> I could use NanoSim, but that tool seems to be very difficult to use. I only have a reference genome, but you need training data as well?

Karel Břinda's fork, NanoSim-H, ships with pretrained models (or "error profiles", as they're called there), eliminating the need for training data.

Minimal usage instructions to run with the default error model (trained on 1D² reads of E. Coli on an R9 flowcell):

```
pip install nanosim-h
nanosim-h your-reference.fasta

```
> 3 votes

---
Tags: nanopore, long-reads, simulated-data
---


# Question
Title: Proteomics data Vs Transcriptomics data?
I want to use either of Proteomics or Transcriptomics data for integrating it into my kinetic model. Before proceeding, I want to know what are the advantages of using either of them so that I could make an informed decision on it!

Many studies have shown that the best we can do is integrating both transcriptomics and proteomics data with our kinetic model, but I've some time constraints and have to proceed with only one of those.

*My effort and findings:* I've found from discussions with researchers that gathering transcriptomics data has an amplification step which increases the chance of finding a particular one whereas gathering proteomics data has no such step but has fragmentation and then rejoining which creates many problems(such as splice variants etc) and thus leads to a loss of data. But a PostDoc told me that even after the loss of that data, I'll get more information from Proteomics data.

I want to know such type of points and want to know if these are valid or not!

# Answer
Ultimately, proteomics data can actually measure the actual concentrations of the enzymes that you are interested, while transcriptomics only measures their rates of change due to ongoing transcription from mRNA.

Rates of change are often treated as proxy for concentrations. This makes sense when you have a constant degradation/dilution environment and are near steady state, since then the rate of change is proportional to the converged concentration. If you're interested in more dynamic environments or believe that modulation of degradation or sequestration might play a significant role in a system, however, that assumption breaks done.

Thus, if you have to choose between proteomics and transcriptomics data, I would suggest proteomics.

As a final note, however, since the link suggests you're interested in metabolites, however, have you consider going directly to metabolomics data? That would let you really get direct information about the particular species of interests.
> 0 votes

# Answer
This is what I found from doing some research. Comments are welcome at any point!

* Capture percentage for data gathering:

+ Transcriptomics Data: There's an amplification step in Transcriptomics data gathering methods. Hence, it's possible to capture almost the totality of the Transcriptome using those methods(scRNA seq methods, Nanopore tech, Spatial transcriptomics, etc)
    + Proteomics Data: There's no amplification step in Proteomics data gathering methods. Those methods have fragmentation and a rejoining step which leads to a loss of the data.
* Scalability of methods:

+ Transcriptomics Data: Scalable methods available
    + Proteomics Data: Less scalable methods for protein studies
* Reference availability:

+ Transcriptomics Data: Fully annotated references available on consortiums such as ensembl biomart, etc
    + Proteomics Data: Universal and Comprehensive Human Proteome reference is still in question and because of the incomplete and inaccurate references, much of the newly generated data is being rendered useless
* Uniformity

+ Transcriptomics Data: Lots of publications on pipelines, and their benchmarking available
    + Proteomics Data: Lack of uniformity across labs/research groups and lack of related literature lead to differences in protein fragmentation and solubilization and differences in algorithms to run analyses
* Technical Bias

+ Transcriptomics Data: Many methods available to tackle cell-bias, noise in data, and batch effects, so as to capture maximum biological variability
    + Proteomics Data: Cell-wide Mass spectrometry has a bias towards the identification of peptides with higher concentration or contamination from other experiments
> 0 votes

---
Tags: transcriptome, networks, modelling, proteomics, multi-omics
---


# Question
Title: How to select only RNA with Hetero atoms from pdb file with python?
I'm trying to separate RNA from protein in a complex protein/RNA PDB file and I want all RNA info with the hetero atoms in between the bases BUT without H20 etc. In short I want RNA part of pdb file without discontinuous lines.

I managed to separate RNA from protein with Bio PDB Select but it consider hetero atoms as amino acid when I use is\_aa(residue). So hetero atoms wont appear in my "only RNA" file.

```
from Bio.PDB import *
from Bio.PDB import PDBParser, PDBIO, Select
import os

class ProtSelect(Select):
    def accept_residue(self, residue):
        return 1 if is_aa(residue) == True else 0

class RNASelect(Select):
    def accept_residue(self, residue):
        return 1 if is_aa(residue) == False and residue.id[0] != "W" else 0

pdb = PDBParser().get_structure("2bh2", "pdb2bh2.ent")
io = PDBIO()
io.set_structure(pdb)
io.save("seqprotest.pdb", ProtSelect())
io.save("seqRNAtest.pdb", RNASelect())

```

EDIT: I managed to do it:

```
class RNASelect(Select):
    def accept_residue(self, residue):
        return 1 if ((residue.id[0] != 'H_ MG') and ((residue.get_resname() == '  A' or residue.get_resname() == '  U' or   residue.get_resname() == '  G' or  residue.get_resname() == '  C' or residue.get_resname() == '  T') or  (residue.id[0][0:2] == 'H_' and (residue.get_resname()[2] == 'U'  or residue.get_resname()[2] == 'A' or residue.get_resname()[2] == 'C' or residue.get_resname()[2] == 'G' or residue.get_resname()[2] == 'T'))))  else 0

```

It's not pretty but it do the job. The idea is to get everything that finish by a nucleotide. I figured out that hetatom's names of RNA chains always finish by a nucleotide letter so i just have to check everything that begins by 'H\_' and ends by 'A,T,C,G,U'. Just remove 'H\_ MG' and you're good.

# Answer
Solution I found: (c/p of body)

```
class RNASelect(Select):
    def accept_residue(self, residue):
        return 1 if ((residue.id[0] != 'H_ MG') and ((residue.get_resname() == '  A' or residue.get_resname() == '  U' or   residue.get_resname() == '  G' or  residue.get_resname() == '  C' or residue.get_resname() == '  T') or  (residue.id[0][0:2] == 'H_' and (residue.get_resname()[2] == 'U'  or residue.get_resname()[2] == 'A' or residue.get_resname()[2] == 'C' or residue.get_resname()[2] == 'G' or residue.get_resname()[2] == 'T'))))  else 0

```

It's not pretty but it do the job. The idea is to get everything that finish by a nucleotide. I figured out that hetatom's names of RNA chains always finish by a nucleotide letter so i just have to check everything that begins by 'H\_' and ends by 'A,T,C,G,U'. Just remove 'H\_ MG' and you're good.
> 3 votes

---
Tags: python, pdb, rna
---


# Question
Title: How to find novel transcripts using GFFcompare?
I am trying to find novel transcripts from an RNA-seq database. Based on the advice I got, it seemed that using Stringtie for transcript assembly is a good way to go, and it supports novel transcript discovery even with the reference GTF file provided (which apparently significantly improves the assembly process). So I tried to follow the protocol provided in the Nature Protocols paper that described the use of Stringtie; the paper also suggested using GFFcompare for comparing the assembled transcripts with the reference, to quantify and find the novel transcripts.

I followed the entire protocol, and this was the output of GFFcompare when comparing the reference GTF with that of one of the samples after mapping -

```
#= Summary for dataset: ./hisat2/ERR188044_chrX.gtf
#-----------------| Sensitivity | Precision  |
        Base level:    51.7     |    79.5    |
        Exon level:    46.7     |    85.2    |
      Intron level:    47.2     |    93.9    |
Intron chain level:    31.2     |    64.4    |
  Transcript level:    31.6     |    52.3    |
       Locus level:    36.6     |    50.1    |

     Matching intron chains:     580
       Matching transcripts:     664
              Matching loci:     397

          Missed exons:    4395/8804    ( 49.9%)
           Novel exons:     426/4874    (  8.7%)
        Missed introns:    3832/7946    ( 48.2%)
         Novel introns:      83/3992    (  2.1%)
           Missed loci:     610/1086    ( 56.2%)
            Novel loci:     273/795 ( 34.3%)

 Total union super-loci across all input datasets: 749
1270 out of 1270 consensus transcripts written in 188044_only.annotated.gtf (0 discarded as redundant)

```

In the paper however, the number of novel **genes** and **transcripts** is clearly mentioned -

So, I was wondering how they were able to calculate the number of novel **transcripts and genes**, from the data that GFFcompare provides (they even mention in the paper that they got it from GFFcompare). I see the numbers related to **novel exons and novel introns**, but how do I translate it to the number and identities of the novel transcripts/genes?

Cross-posted in Biostars. Let me know if I shouldn't do this.

# Answer
> 4 votes
The output of gffcompare includes several files per run (just like cuffcompare). Example for a run:

```
$ ls cuffcmp.* | sed 's/\t/\n/'
cuffcmp.combined.gtf
cuffcmp.loci
cuffcmp.output.gtf.refmap
cuffcmp.output.gtf.tmap
cuffcmp.stats
cuffcmp.tracking

```

From my experience, the easiest file to manipulate is the `tmap` file:

> Tab delimited file lists the most closely matching reference transcript for each query transcript. There is one row per query transcript, and the columns are as follows:
>
> ```
> 1     Reference gene name    The gene_name attribute of the reference GTF record for this transcript, if present. Otherwise gene_id is used.
> 2     Reference transcript id     The transcript_id attribute of the reference GTF record for this transcript
> 3     Class code  The type of relationship between the query transcripts in column 4 and the reference transcript (as described in the Class Codes section below)
> 4     Query gene id   The query (e.g., Stringtie) internal gene id
> 5     Query transcript id     The query internal transcript id
> 6     Number of exons     The number of exons in the query transcript
> 7     FPKM    The expression of this transcript expressed in FPKM
> 8     TPM     the estimated TPM for the transcript, if found in the query input file
> 9     Coverage    The estimated average depth of read coverage across the transcript.
> 10    Length  The length of the transcript
> 11    Major isoform ID    The query ID of the gene's major isoform
> 12    Reference match length  The length of the longest overlap with a reference, ‘-’ if there is no such exonic overlap
>
> ```

The important column to look at for transcript classification is actually the 3rd, "class code". The content is meaningful if you used option `-r` when running gffcompare.

In the tmap file you can access info on both transcript and gene ID directly, both from the assembler and the reference annotation.

Now what you need is:

1. n. assembled genes
2. n. novel genes
3. n. novel transcripts
4. n. transcripts matching annotation

At this point, it is essential to understand what the authors mean by "novel" (gene/transcript), does it mean "non previously annotated" (= class code `u`, for "unknown")? Also, what about "matching", does it mean "perfectly matching" (= class code `=`, for "complete, exact match of intron chain")?

Under these assumptions, you can do some easy stats with some bash help.

1. n. assembled genes:

`$cat cuffcmp.output.gtf.tmap | sed "1d" | cut -f4 | sort | uniq | wc -l`
2. n. novel genes:

`$ cat cuffcmp.output.gtf.tmap | awk '$3=="u"{print $0}' | cut -f4 | sort | uniq | wc -l`
3. n. novel transcripts:

`$ cat cuffcmp.output.gtf.tmap | awk '$3=="u"{print $0}' | cut -f5 | sort | uniq | wc -l`
4. n. transcripts matching annotation:

`$ cat cuffcmp.output.gtf.tmap | awk '$3=="="{print $0}' | cut -f5 | sort | uniq | wc -l`

I hope this helps. I think it's an easier approach than trying to "translate the numbers related to novel exons and novel introns to the number and identities of the novel transcripts/genes" from the stats file.

Although you have to think a bit what "class codes" to include in each count. Maybe my assumptions are too strict for your use case.

---
Tags: rna-seq, transcriptome, reference-genome
---


# Question
Title: DoHeatmap Hierarchical Clustering Seurat
I have a list of genes that I'd like to visualize using the DoHeatmap function in Seurat. However, the output of the heatmap does not result in hierarchical clustering and therefore makes it very difficult to interpret. My code is,

```
 DoHeatmap(object = obj, genes.use = genes),
          slim.col.label = TRUE, remove.key = T,group.label.rot = F, use.scaled = T)

```

Is there a way to adjust the `DoHeatmap` command to rank the cells by the intensity of gene expression? Here's an example output:

# Answer
> 1 votes
You need to order the marker matrix (e.g. by `avg_logFC`) before calling `DoHeatMap`.

```
library(dplyr)
all.markers <- FindAllMarkers(object = obj)
top20 <- all.markers %>% group_by(cluster) %>% top_n(20, avg_logFC)
DoHeatmap(object = obj, genes.use = top20$gene, slim.col.label = TRUE, remove.key = TRUE)

```

---
Tags: r, seurat, heatmap
---


# Question
Title: Alternative to enrichR for enrichment analysis?
I have been using EdgeR from the **Bioconductor** package to perform differential expression analysis on bulk RNASeq data sets.

This analysis gives me an output dataframe that has a list of all the genes, ranked by the their adjP values. This dataframe also contains the t-statistic and lfc values for each gene as well.

I have recently been using enrichR to perform gene set enrichment analysis on all the differentially expressed genes to determine what biological, molecular, transcriptional, etc pathways these genes are involved in. I really like the enrichR package for a few key reasons:

1.) It has access to many different databases that I can call on with a simple vector of their names. Most of my data is murine RNASEq, so for those analyses, I can simply select the murine databases:

```
enrichR_dbs <- c("WikiPathways_2019_Mouse", "KEGG_2019_Mouse", "GO_Biological_Process_2018", "GO_Cellular_Component_2018", "GO_Molecular_Function_2018", "GWAS_Catalog_2019", "TRRUST_Transcription_Factors_2019")

```

2.) It uses a single function, enrichR, to perform the enrichment. The arguments for that function are the list of gene symbols and the list of databases.

```
enrichr(as.character(Gene_Dataframe$SYMBOL[1:500]), enrichR_dbs)

```

Anyway, I am looking for another enrichment analysis because I want to be able to see if I get similar pathways for my genes when using another method/algorithm.

Unfortunately, I can't seem to find another enrichment analysis( competitive or self-contained) that is both as user-friendly and that gives me easy access to murine databases (without requiring to convert my genes through changes in capitalization, like you need to with GSEA) as enrichR.

**Does anyone have recommendations about enrichment analysis R packages that don't require too many more inputs than the dataframe I generated through my EdgeR analysis?**

P.S. I'm ok with downloading databases if I need to. Though inclusion of a PDF with a tutorial would be fantastic, if possible.

# Answer
A couple of years ago I wrote a blog post about GSEA methods in Bioconductor, in addition to these I would recommend clusterProfiler or ReactomePA. In most of them you need to provide the gene set/pathway database.
If you want to do gene ontology enrichment then I recomment the topGO package. In most, if not all, of these methods you just need an ordered list and a gene set source.
> 2 votes

# Answer
Try gProfileR. It is very intuitive and easy to use in R.
> 2 votes

---
Tags: go-enrichment
---