text:
How can I assign a genomic region into a window using R?
I have asked this question at BioStars but have not gotten any suggestions so far. So, I am asking the community here, perhaps someone here has an idea how I can solve my problem.
I am dealing with a bunch of inversions and running some analysis on each of them. One of them is located on chromosome one, the starting position is 6199548 and ending position is 9699260. So the length of this variant is 3499712 bp. I want to generate non-overlapping windows based on the length of these focal positions (3499712) and assign the focal position to one of these windows that overlap with the starting and ending position of my focal position. I am achieving this by using the GenomicRanges::tileGenome function in R:
seqlength <-c(chr1=159217232,chr2=184765313,chr3=182121094,chr4=221926752,chr5=187413619)
bin_length<-c(3499712)
 bins<-tileGenome(seqlength,tilewidth=bin_length,cut.last.tile.in.chrom
        = T)
            #bins_u<-unlist(bins)
            write.table(bins,file=paste("bins_good_",bin_length,sep=""),col.names
        = FALSE, row.names = TRUE) }

This generates non-overlapping 3499712 bps windows across the genome.
window_id chromosome start end length
"1" "chr1" 1 3499712 3499712
"2" "chr1" 3499713 6999424 3499712
"3" "chr1" 6999425 10499136 3499712
"4" "chr1" 10499137 13998848 3499712
"5" "chr1" 13998849 17498560 3499712

However, I have a major problem here that I have not been able to sort out. I want one of these windows to overlaps exactly with the starting and ending position of my focal position (6199548:9699260). Only the focal position not any other position outside this region. This means 2 windows (one before and one after the focal position) will have a different length which my analysis is fine with. I just do not know how to do it?
This is my desired output:
window_id chromosome start end length
"1" "chr1" 1 3499712 3499712 "*"
"2" "chr1" 3499713 6199547 2699834 "*"
"3" "chr1" 6199548  9699260 3499712 "*"
"4" "chr1" 9699261 13998848 4299587 "*"
"5" "chr1" 13998849 17498560 3499712 "*"

As you see, the the second window is now smaller, the third window accommodates the focal position and the forth window is a little larger than the bin size. Does anyone have any idea how can I do this task in R?

substr(chromosome, start pos, finish pos)
?

I'm not sure how you did your calculation but for me it looks like your region is 3499713 bp in length.
site <- GRanges("chr1:6199548-9699260")

seqlength <- c(chr1=159217232,chr2=184765313,chr3=182121094,chr4=221926752,chr5=187413619)

bin_length <- width(site)

## check length of region
print(bin_length)
##3499713

bins <- tileGenome(seqlength, tilewidth=bin_length, cut.last.tile.in.chrom=TRUE)

## extract overlapping bins
flank.bins <- subsetByOverlaps(bins,site)

## find new coordinates excluding region of interest
new.flank.bins <- setdiff(reduce(flank.bins),site)

## combine new bins with your region of interest
new.bins <- c(new.flank.bins,site)

## remove old bins
bins <- subsetByOverlaps(bins,new.bins,invert=T)

## add new bins
bins <- sort(c(bins,new.bins))

## check that it makes sense
##table(width(bins[seqnames(bins) == 'chr1']))

write.table(bins,file=paste("bins_good_",bin_length,sep=""),col.names= FALSE, row.names = TRUE)

I sincerely appreciate your help. I did not know that the GenomicRanges options! I am now working on it to learn about the many other useful functions.

I have been digging into your code, but I guess it does not give the outputs that I want. This is the output that I get from your code:"1" "chr1" 1 3499713 3499713 ""
"2" "chr1" 3499714 10499139 6999426 ""
"3" "chr1" 6199548 9699260 3499713 ""
"4" "chr1" 10499140 13998852 3499713 "". this is what I get , the focal window(window 3) overlaps with window 2!

If you copy and paste my code exactly and look at the first 5 lines of the resulting file, you should get the following: "1" "chr1" 1 3499713 3499713 ""
"2" "chr1" 3499714 6199547 2699834 ""
"3" "chr1" 6199548 9699260 3499713 ""
"4" "chr1" 9699261 10499139 799879 ""
"5" "chr1" 10499140 13998852 3499713 "*" If you share your code I can help you further.


text:
BiomaRt error: Error in martCheck
I am experiencing an error while using getSequence() in biomaRt.
The code is given below:
library("biomaRt")
ensemblp <- useMart(biomart="plants_mart",host="plants.ensembl.org")
osj = useDataset("osativa_eg_gene", mart = ensemblp)
genes <- read.csv("osj.csv")
osj_id <- genes$ensembl_transcript_id

os_genes <- getSequence(id = osj_id, type = "ensembl_transcript_id", seqType = "peptide", mart = ensemblp)

The error message I am getting is :
Error in martCheck(mart, c("ensembl", "ENSEMBL_MART_ENSEMBL")) :
This function only works when used with the ensembl BioMart.This
function only works when used with the ENSEMBL_MART_ENSEMBL BioMart.

Did I make any mistake in code? Please help me fix the error. I tried asking this in Bioconductor support, but no one responded yet.

It doesn't seem to be a mistake in the code, but simply what you try to do isn't working in that BioMart

Did you find and review this post on bioconductor?  Does it help you? https://support.bioconductor.org/p/74761/

getSequence has only been enabled for the main Ensembl (vertebrates) biomaRt. This is all at the Bioconductor end. You'll need to use a getBM instead, for example:
os_genes <- getBM(attributes =c("ensembl_transcript_id", "peptide"), filters = "ensembl_transcript_id", values = osj_id, mart = osj)

Biomart seems to be down most of the time so i think you can parse the ensemble ID and names from the annotation file which i did as i was having similar problem .
You can use this to make a file with ensemble ID and gene names for you use
awk '{
    for (i = 1; i <= NF; i++) {
        if ($i ~ /gene_id|gene_name/) {
            printf "%s ", $(i+1)
        }
    }
    print ""
}' Homo_sapiens.GRCh37.70.gtf | sed -e 's/"//g' -e 's/;//g' -e 's/ /\t/' | sort -k1,1 | uniq > Homo_sapiens.GRCh37.70.txt

source: Stack Exchange

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

sort: votes


text:
How to output all sequences with bwa mem, not `*`?
I've been running bwa mem -a for alignment, using the -a flag---this will

output all alignments for SE or unpaired PE

I've noticed in the SAM that there are several alignments with * in the SEQ and QUAL fields. Based on the documentation:

SEQ: segment SEQuence. This field can be a ‘*’ when the sequence is not stored. If not a ‘*’, the length of the sequence must equal the
  sum of lengths of M/I/S/=/X operations in CIGAR. An ‘=’ denotes the
  base is identical to the reference base. No assumptions can be made on
  the letter cases.
QUAL: ASCII of base QUALity plus 33 (same as the quality string in the Sanger FASTQ format). A base quality is the phred-scaled base
  error probability which equals −10 log10 Pr{base is wrong}. This field
  can be a ‘*’ when quality is not stored. If not a ‘*’, SEQ must not be
  a ‘*’ and the length of the quality string ought to equal the length
  of SEQ.

it appears the sequence isn't stored.
I would strongly prefer to have the sequence in this case. Is there any to direct bwa to output these sequences?

Are you aware that many of these records are secondary/supplementary alignments? Does your workfow/model require/benefit from these alignments and handle them correctly?

@DanielStandage "Are you aware that many of these records are secondary/supplementary alignments?" Yes. " Does your workfow/model require/benefit from these alignments and handle them correctly?". Benefit, yes.

I don't believe that is an option.  Internally we created a patch and submitted it to BWA however that was 4 years ago and it has not been accepted.

Thanks. Did @user172818 ever give a reason why this feature isn't included?

see also https://www.biostars.org/p/352930/

@Pierre This is a great tool Pierre! I think it runs successfully for me. The problem is, I get errors downstream when using samtools view, i.e.. [E::sam_parse1] CIGAR and query sequence are of different length Are there work arounds?

I didn't test it with all the possible cases. Can you please try to isolate a minimal sam file containing the faulty reads ( all the reads with the same name) and send a bug report ? https://github.com/lindenb/jvarkit/issues/

@Pierre I may send you an e-mail instead. Thanks Pierre!

plindenbaum reverse(oohay) fr but I'll need a minimal failing SAM file.

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


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

Do you really mean bedpe? I may be missing something, but what would converting a vcf to bedpe even mean? Bedpe is normally used to described paired sequences, usually primers. What would converting a VCF to bedpe entail? Could you maybe [edit] your question and give an example of a minimal VCF and what it would look like as bedpe?

@terdon I've edited the question for clarity

Thanks! That makes it much clearer. I don't really see how you'd represent deletions, but if the LUMPY folks are doing it, I guess it must be possible. Did you try vftobedpe from svtools?

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

Could you edit your answer and add this command's name and how to execute it? As it stands, this is kind of more of a comment than an answer.

Feel free to move it to a comment then. I consider the documentation of the tool sufficient and clear, and do not believe that copy-pasting that usage information is beneficial, as it might get outdated while the official documentation shouldn't.

It's a question of what is considered an answer on the SE platform. This seems like a good answer, but we really want the relevant information to be here in case the link goes stale and to be able to have the relevant details here. You're linking to the docs, so if they change and your instructions become obsolete, the docs will still be there. In the meantime, however, this is not actually giving us enough info to know how to use it and that's a shame.

Hi @terdon, thanks for the edits. "However, this isn't a bedpe file, so it looks like it isn't what the OP was asking for". That looks pretty close to a BEDPE, right? What is the difference you see, besides the file extension?

Whoops, my bad! It is indeed a bedpe and not a bed. I had managed to miss the second set of coordinates, sorry! However, it doesn't seem to be right for what the OP wants since that's just repeating the start position two times and then the end position another two times. Have a look at the blog post the OP linked to for what they seem to be looking for.

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


text:
What does an output of -I for all amino acids mean in a psi-blast pssm matrix file?
I have run psi-blast using the NR database, remotely, with one iteration, for several sequences to calculate an evolutionary profile (PSSM) for each of those sequences.
However, many of the PSSM files contain lines of -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I rather than numerical values.  What does it mean?
>     Last position-specific scoring matrix computed
>                 A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V
>         1 G  -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I
>         2 N     4  -1  -2  -2   0  -1  -1   0  -2  -1  -1  -1  -1  -2  -1   1   0  -3  -2   0
>         3 K    -2  -1   4   4  -3   0   1  -1   0  -3  -4   0  -3  -3  -2   0  -1  -4  -3  -3
>         4 E     0  -3  -3  -3   9  -3  -4  -3  -3  -1  -1  -3  -1  -2  -3  -1  -1  -2  -2  -1
>         5 K    -2  -2   1   6  -3   0   2  -1  -1  -3  -4  -1  -3  -3  -1   0  -1  -4  -3  -3
>         6 A    -1   0   0   2  -4   2   5  -2   0  -3  -3   1  -2  -3  -1   0  -1  -3  -2  -2
>         7 D    -2  -3  -3  -3  -2  -3  -3  -3  -1   0   0  -3   0   6  -4  -2  -2   1   3  -1
>         8 R     0  -2   0  -1  -3  -2  -2   6  -2  -4  -4  -2  -3  -3  -2   0  -2  -2  -3  -3
>         9 Q    -2   0   1  -1  -3   0   0  -2   8  -3  -3  -1  -2  -1  -2  -1  -2  -2   2  -3
>        10 K    -1  -3  -3  -3  -1  -3  -3  -4  -3   4   2  -3   1   0  -3  -2  -1  -3  -1   3
>        11 V    -1   2   0  -1  -3   1   1  -2  -1  -3  -2   5  -1  -3  -1   0  -1  -3  -2  -2
>        12 V    -1  -2  -3  -4  -1  -2  -3  -4  -3   2   4  -2   2   0  -3  -2  -1  -2  -1   1
>        13 S    -1  -1  -2  -3  -1   0  -2  -3  -2   1   2  -1   5   0  -2  -1  -1  -1  -1   1
>        14 D    -2   0   6   1  -3   0   0   0   1  -3  -3   0  -2  -3  -2   1   0  -4  -2  -3
>        15 L    -1  -2  -2  -1  -3  -1  -1  -2  -2  -3  -3  -1  -2  -4   7  -1  -1  -4  -3  -2
>        16 V    -1   1   0   0  -3   5   2  -2   0  -3  -2   1   0  -3  -1   0  -1  -2  -1  -2
>        17 A    -1   5   0  -2  -3   1   0  -2   0  -3  -2   2  -1  -3  -2  -1  -1  -3  -2  -3
>        18 L     1  -1   1   0  -1   0   0   0  -1  -2  -2   0  -1  -2  -1   4   1  -3  -2  -2
>        19 E     0  -1   0  -1  -1  -1  -1  -2  -2  -1  -1  -1  -1  -2  -1   1   5  -2  -2   0
>        20 G     0  -3  -3  -3  -1  -2  -2  -3  -3   3   1  -2   1  -1  -2  -2   0  -3  -1   4
>        21 A    -3  -3  -4  -4  -2  -2  -3  -2  -2  -3  -2  -3  -1   1  -4  -3  -2  11   2  -3
>        22 L    -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
>        23 D    -2  -2  -2  -3  -2  -1  -2  -3   2  -1  -1  -2  -1   3  -3  -2  -2   2   7  -1
>        24 M    -1   0   0   1  -3   4   4  -2   0  -3  -3   1  -1  -3  -1   0  -1  -2  -2  -2
>        25 Y     0  -3  -3  -3   9  -3  -4  -3  -3  -1  -1  -3  -1  -2  -3  -1  -1  -2  -2  -1
>        26 K    -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4
>        27 L    -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
>        28 D    -1  -2  -3  -3  -1  -2  -3  -4  -3   3   3  -3   2   0  -3  -2  -1  -2  -1   2
>        29 N  -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I
>        30 S    -2   1  -1  -1  -2  -1  -2  -2  -1  -1  -1  -1  -2   0  -1   0   0  -1  -1  -3
>        31 R     4   0  -2   1   4  -1  -3   0  -3   0  -3  -4   4  -1   0  -1  -3  -1   0  -4
>        32 Y    -3  -1  -3  -4  -3  -3  -2  -3  -1  -3  -1  -1  -3  -3  -3  -1  -1  -1  -3  -2

I don't know but guess that either something has gone wrong or those somehow represent insertions/deletions (gaps)

can you try running with more iterations? Running with a single one makes little sense I think

@Chris_Rands Perhaps more iterations would help, however, remote blast only allows one iteration, unless requerying.  I could try running locally to see.  Strangely, for SWISSPROT database (tiny compared to NR) it doesn't appear to happen (the instances of -I).

It means -inf. It has to do with how the PSSM score is calculated.

That sounds very plausible, do you have a source/reference? Aren't pseudocounts normally added by default? will running >1 iteration change this?

Yes it appears to be minus infinity according to that wikipedia article.  Strange how it was so hard to find this information.  Thanks!  Now I just need to figure out what to do about it since I need to encode this data.

@Chris_Rands, if you have blast installed type in psiblast -help. You'll see that pseudocount is set by default to 0.

On NCBI PSI-BLAST web interface it says "Pseduocount parameter. If zero is specified, then the parameter is automatically determined through a minimum length description principle (PMID 19088134). A value of 30 is suggested in order to obtain the approximate behavior before the minimum length principle was implemented."

@Chris_Rands, did you read psiblast -help or not? That is something else than a web interface.


text:
How can I only get the species name for fasta sequences from blast results?
I am trying to make a phylogenetic tree from sequences obtained with blast. I have files that contains hundreds of fasta protein sequences which are all named like >NP_567483.1 transcriptional regulator [Pseudomonas aeruginosa PAO1]. But if I use these files to do multiple alignments or phylogenetic trees using either ClustalX or MEGA7, every sequence are named after the first ID in the sequence names (for this example it would be NP_567483.1), and what I want is to have the name of the strain instead (here Pseudomonas aeruginosa PAO1) which is part of the sequence description.
So I don't know whether the softwares I'm using understand that the information between the [ ] is the species and that I just didn't find where to set this as the name or if I have to write some kind of a script to change my fasta files before. I'm really new to this and don't know much about it. Thanks for your help in advance. Julian

Should we assume one sequence per species? Otherwise species name is no good as an identifier.

There might be a "cleaner" way to do this, but what I would do is first remove all of the spaces from the fasta sequence names and replace them with underscores using something like sed.
sed 's/ /_/g'

Run your BLAST after you remove the spaces and you will get the whole ID string in the result instead of only that first part (e.g. NP_567483.1). From here you can parse out the species name, or whatever other ID you would like.
Hope this is helpful! :)

You can replace fasta header of each sequence with the string that is inside the brackets using awk:
awk -F '[][]' '{if ($0~/>/) {print ">"$2} else {print $0}}' input.fasta > output.fasta

For lines starting with > (fasta header), it replaces the line with the string inside the brackets and for other lines (which contain the sequence), it prints the line without making any changes.

With sed, first test it:
echo '>NP_567483.1 transcriptional regulator [Pseudomonas aeruginosa PAO1]' | sed "s/>.*\[\(.*\)\]/>\1/"

Then modify your fasta:
sed -i.bak "s/>.*\[\(.*\)\]/>\1/" file.fasta

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


text:
HMM Profile from Convergence of retrotransposons in oomycetes and plants
I am trying to recreate the experiment from Convergence of retrotransposons in
oomycetes and plants by Kirill Ustyantsev, Alexandr Blinov and Georgy Smyshlyaev. In the "Methods" section a

HMM profile was constructed from the amino acid alignment of aRNH
  sequences from the Ustyantsev et al

I was wondering where I could get this profile. I checked the cited papers but they were also referencing other papers. Could someone maybe clarify where I could get the profile?

Your question is how to to this HMM or how to obtain/download it? In the later case you can always try to contact the authors if it is important to use exactly the same HMM profile

The alignment they used to build the HMM is included as supplementary file 8 in the paper you cite:

Multiple amino acid sequence alignment of RNH genes and domains from diverse taxa constructed and used for the phylogenetic reconstruction in the present study. (TXT 45 kb)

So, you can download that alignment and then use a tool like hmmbuild (part of HMMER) to build the HMM profile:
hmmbuild model.hmm 13100_2017_87_MOESM8_ESM.txt

The command above will create the file model.hmm which contains the HMM built from the alignment in the supplementary file 13100_2017_87_MOESM8_ESM.txt. I am including the first few lines here (this is not the entire file, don't use this as the model—build it yourself. I only include it to illustrate the format):
HMMER3/f [3.1b2 | February 2015]
NAME  13100_2017_87_MOESM8_ESM
LENG  121
ALPH  amino
RF    no
MM    no
CONS  yes
CS    no
MAP   yes
DATE  Mon Dec 11 13:33:29 2017
NSEQ  165
EFFN  5.697556
CKSUM 2893352190
STATS LOCAL MSV       -9.8975  0.71332
STATS LOCAL VITERBI  -10.6450  0.71332
STATS LOCAL FORWARD   -3.9247  0.71332
HMM          A        C        D        E        F        G        H        I        K        L        M        N        P        Q        R        S        T        V        W        Y
            m->m     m->i     m->d     i->m     i->i     d->m     d->d
  COMPO   2.34473  4.21240  2.85435  2.69203  3.54228  2.79153  3.54812  3.01140  2.72630  2.48370  3.82225  2.94269  3.49156  3.05882  2.84632  2.63343  2.86306  2.72633  4.27012  3.56181
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.13773  5.45955  2.08419  0.61958  0.77255  0.00000        *
      1   2.99488  4.57852  4.44659  3.29307  1.80597  4.16203  3.88788  2.87265  3.75733  2.38594  3.30344  3.33633  4.53284  3.64324  3.69630  2.90350  3.30742  2.83413  4.36991  1.50673      1 y - - -
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.00723  5.32905  6.05140  0.61958  0.77255  0.45275  1.01027
      2   2.66602  2.17750  5.01715  4.40788  1.72617  4.29647  4.62737  2.89579  4.21116  2.71668  3.64495  4.43195  4.66188  4.35829  3.99238  3.10363  2.02207  1.89528  5.14692  3.02155      2 f - - -
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.00661  5.41881  6.14115  0.61958  0.77255  0.61842  0.77391

thanks for the clarification. I missed the additional file since I was looking at a printed out copy.


text:
What is the meaning of the Hmmer profile?
I want to understand what the hmmer Profile used in hmmer.
consider the profile
HMMER2.0  [converted from 3.1b2]
NAME  aRNH_profile
LENG  121
ALPH  Amino
RF    no
CS    no
MAP   yes
NSEQ  34
DATE  Tue Apr 26 21:24:57 2016
XT       -8455     -4  -1000  -1000  -8455     -4  -8455     -4
NULT        -4  -8455
NULE       656  -1722     98    419   -333    475  -1125    239    250    947  -1073   -271    -50   -338    114    451    113    430  -2131   -717
HMM        A      C      D      E      F      G      H      I      K      L      M      N      P      Q      R      S      T      V      W      Y
         m->m   m->i   m->d   i->m   i->i   d->m   d->d   b->m   m->e
           -6      *  -7932
     1    -214   -916   -565    694    862  -1448    512   -838    115  -1042    170   1348  -1481    369   -158   -290   -159   -625    -52   1675     1
     -    -210   -336    220    -37   -342    377     73   -669    210   -511   -732    404    423     70     27    439    205   -414   -161   -176
     -     -18  -6890  -7932   -894  -1115   -701  -1378    -18      *
     2     285    990  -2213  -1678   2305  -1819   -698    219  -1289   -140    545  -1167  -1835   -973  -1322   -424   1025    762   -559    377     2
     -    -210   -336    220    -37   -342    377     73   -669    210   -511   -732    404    423     70     27    439    205   -414   -161   -176
     -     -18  -6890  -7932   -894  -1115   -701  -1378      *      *
[....]

what exactly does each component mean? By definition a HMM $\lambda$ consist of an alphabet $V$ in my case twenty amino acids. Different states $S$ in our case main, insertion and deletion (beginning and end). And two probability matrices of state transitions and emission in certain state.
now if I look at the profile I can see that it has the allowed transitions
m->m   m->i   m->d   i->m   i->i   d->m   d->d   b->m   m->e
and all the amino acids, but what does each row in the table mean. I checked the manual from the documentation and it says rather ironically

If you were to look at this file (and you don’t have to – it’s
  intended for HMMER’s consumption, not yours), you’d see [part of my output]

are these number translations of the probabilities in $[0,1]$  into another interval $[-x,+y]$? what is the meaning of the three rows under A 1
-124
-210
 -18

Is this the transition matrix? It would be great if someone could explain this  or point to the correct reference.

The tutorial in the documentation does indeed state what you quoted but a full description of the file format is given in Section 8 of the documentation, starting on page 106.
To summarise, as you noted, the number columns A-K is the alphabet size (the number of amino acids). The number of rows underneath is the number of nodes, which corresponds to the number of match states. The three numbers (-214, -210 and -18 in your example) are: the match emission, insert emission, and state transition values.
Caveat: The HMMER profile HMM format changes frequently between HMMER versions (even between minor versions) so there is no guarantee this information is correct for all the older file formats.

Quoting the relevant part of the documentation on page 109:

main model section
All the remaining fields are mandatory. The first
two lines in the main model section are atypical.‡ They contain
information for the core model’s BEGIN node. This is stored as model
node 0, and match state 0 is treated as the BEGIN state. The begin
state is mute, so there are no match emission probabilities. The first
line is the insert 0 emissions. The second line contains the
transitions from the begin state and insert state 0. These seven
numbers are: B → M1, B → I0, B → D1; I0 → M1, I0 → I0; then a 0.0 and
a ’’, because by convention, nonexistent transitions from the
nonexistent delete state 0 are set to log 1 = 0 and log 0 = −∞ = ‘’.
The remainder of the model has three lines per node, for M nodes
(where M is the number of match states, as given by the LENG line).
These three lines are (K is the alphabet size in residues):
Match emission line
The first field is the node number (1 . . . M). The
parser verifies this number as a consistency check (it expects the
nodes to come in order). The next K numbers for match emissions, one
per symbol, in alphabetic order. The next field is the MAP annotation
for this node. If MAP was yes in the header, then this is an integer,
representing the alignment column index for this match state
(1..alen); otherwise, this field is ‘-’. The next field is the CONS
consensus residue for this node. If CONS was yes in the header, then
this is a single character, representing the consensus residue
annotation for this match state; otherwise, this field is ‘-’. The
next field is the RF annotation for this node. If RF was yes in the
header, then this is a single character, representing the reference
annotation for this match state; otherwise, this field is ‘-’. The
next field is the MM mask value for this node. If MM was yes in the
header, then this is a single ’m’ character, indicating that the
position was identified as a masked position during model
construction; otherwise, this field is ‘-’. The next field is the CS
annotation for this node. If CS was yes, then this is a single
character, representing the consensus structure at this match state;
otherwise this field is ‘-’.
Insert emission line
The K fields on this
line are the insert emission scores, one per symbol, in alphabetic
order.
State transition line
The seven fields on this line are the
transitions for node k, in the order shown by the transition header
line: Mk → Mk+1, Ik, Dk+1; Ik → Mk+1, Ik; Dk → Mk+1, Dk+1. For
transitions from the final node M, match state M + 1 is interpreted as
the END state E, and there is no delete state M + 1; therefore the
final Mk → Dk+1 and Dk → Dk+1 transitions are always * (zero
probability), and the final Dk → Mk+1 transition is always 0.0
(probability 1.0).

Thanks, Chris. I missed the later chapters. Sorry that was too quick of me.


text:
SingleCellExperiment: invalid rownames length
I did all of the installations from scratch in a virtual conda environment. Then I tried to run R script from within javascript in the following way:
let shCMD = 'Rscript --vanilla R/full_analysis.R ';
  shCMD += '--diff_1 ' + diff_str_1 + '--diff_2 ' + diff_str_2;
  shCMD += '--markers ' + filterStr + '--markers-highlighted ' + colorStr;
  shCMD += '--datasets ' + datasetStr;

  console.log('constructed command for full analysis: ');
  console.log(shCMD);
  shell.exec('source activate rmain', function(env_code, env_stdout,
env_stderr) {
    console.log('activated rmain environment')
    if (env_code === 0) {
      shell.exec(shCMD, function(code, stdout, stderr) {

As can be seen I am activating the environment by source activate rmain and I am pretty sure that the next shell command that invokes R script runs in that environment because R script needs argparse package and it is installed only in the rmain virtual environment: in rmain argparse package is successfully imported, but outside of it it gives an error for the argparse absence.
So, now the following line in the R script:
sce <- SingleCellExperiment(list(counts=as.matrix(data[,-1])))

gives me the error:
Error in `rownames<-`(`*tmp*`, value = c("0610007P14Rik", "0610009B22Rik",  :
  invalid rownames length
Calls: SingleCellExperiment ... elementMetadata -> .local -> rownames<- -> rownames<-
Execution halted

However, if I launch Rstudio from within the rmain environment, the line runs successfully... So, I am puzzled now on what is going on.

If it works interactively (in Rstudio) but not as a script it could be because R loads certain libraries for you in an interactive session. You could try adding library(methods) to the script and see if you still have the same problem


text:
Question regarding the function pairwiseAlignment in R
I'm a bit confused and would appreciate your assistance, please.
I'm trying to use the pairwiseAlignment function from the Biostrings package to align two protein sequences.
I'm not quite sure where to put my reference sequence, the sequence I want to compare to. Is it under "subject" or "pattern"?
Also, is the score I get the edit distance? Is it the alignment score? Or is it the the alignment length? When I use local alignment, will the score be higher the more the sequences resemble each other? And when I use global alignment, will the score be lower the more the sequences resemble each other?

For semantic purposes, it's better to use the reference sequence as the pattern. Note that only the pattern argument can be used as a multi-element vector, i.e. you can match multiple patterns against a single subject with this function, but not a single pattern against multiple subjects. Knowing that might alter what is defined as a "pattern" or a "subject".
The scores are associated with the alignment method that is used. Alignment scores are higher for more similar matches and lower for poor matches, but the weightings associated with matches and mismatches dictate what that actually means. It might be that a score of -2 happens to be the optimal alignment score, because all other possible alignments are more negative.
It's possible to use pairwiseAlignment to generate the levenshtein distance (most commonly-used edit distance) between two sequences by setting the match score to 0, the mismatch score to -1, gap opening penalty to 0, gap extension penalty to 1, and then negating the resultant pairwiseAlignment score.

According to the help page, it seems more logic to put the reference sequence in subject, and your query sequence in pattern, but it won't really matter for the outcome (except that the pattern sequence will be on top of the alignment).
For your second question:

the scores for the optimal pairwise alignments is returned


text:
Identifying relevant SNPs from a list
I have a list of all SNPs (a list of RsIDs) that fall in the coding sequence of several thousand human genes, and I'm looking for a convenient method to come up with a list of known SNPs that alter an amino acid in those genes.
I'll be glad to hear your thoughts.

Are you open to commercial (non-free) solutions? Also, http://varso.me is a great tool but it can't take batch requests.

What you are looking for is SNP annotation. If you have the chromosome:position reference and alternate alleles for your SNPs of interest, it can be as simple as uploading them to the variant effect predictor.
http://grch37.ensembl.org/Homo_sapiens/Tools/VEP
This will give you the predicted protein change and novelty of the variant with respect to known variation in large scale sequencing projects like the 1000 genomes.
If you cannot afford to upload data on the web, there are downloadable command line software like SNPeff, Variant effect predictor etc but it would require some more bioinformatics knowledge of converting your data to a VCF format and obtaining a human reference genome.

ANNOVAR is another tool that will help you functionally annotate genetic variants. It will tell you if the variant is known in COSMIC or Clinvar databases, which AA change will occur, if it's a frameshift, synonymous, non-synonymous, etc...


text:
UCSC liftover cannot produce bed format output
I was trying to use UCSC liftover tool to convert human coordinates to mouse. While most of the files worked, one file didn't:
chr16   86430087    86430726    Human_chr16:86430087-86430726___element_1___positive____neural_tube[12/12]___hindbrain_(rhombencephalon)[12/12]___limb[3/12]___cranial_nerve[8/12]
chr16   85620095    85621736    Human_chr16:85620095-85621736___element_2___negative_
chr16   80423343    80424652    Human_chr16:80423343-80424652___element_3___negative_
chr16   80372593    80373755    Human_chr16:80372593-80373755___element_4___positive____neural_tube[6/10]___hindbrain_(rhombencephalon)[10/10]___midbrain_(mesencephalon)[10/10]

Does anyone know why this bed file can only be converted to chr:start-end format but not bed4 format?

Were there any logged messages during the process?  Can you include an example that did complete the liftover?

If you remove everything after the "Human" in each line then it works

I figured out the issue by myself. It turned out that UCSC liftover doesn't like the character ":"
I removed ":" and it worked


text:
How to filter ribosomal RNA from scRNA-seq data
I want to filter out ribosomal RNA from scRNA-seq data (downloaded from here). Is there a list of known ribosomal RNA?
The only solution I found is SortMeRNA, however it works with raw sequencing data afaik, while I already have a matrix with transcript counts for each cell. I searched for a comprehensive list of rRNAs but I didn't find any.

What GTF file did you use to generate the counts? For some species only 5S (and maybe MT) rRNA is left, for others you might have one or more copies of the various 48S components.

I actually don't have GTF files. I downloaded the data provided by the Tabula Muris Consortium (https://figshare.com/projects/Tabula_Muris_Transcriptomic_characterization_of_20_organs_and_tissues_from_Mus_musculus_at_single_cell_resolution/27733)

The rRNA genes in that dataset are Rn45s and Rn4.5s.
BTW, you have gene counts, not transcript counts.

Thanks. What do you mean by gene counts instead of transcript counts? That the counts are already collapsed by gene?

Also, how did you find the list of rRNA genes in such dataset?

@gc5 The counts were likely only ever done by genes to begin with. rRNA genes typically have names starting with "Rn" in mouse.

Ok. Should I need to find ribosomal genes in human do you know any database which I can use?

I think they have similar names in humans.

In the paper mentioned, we used the ScaleData function in Seurat to regress out the number of reads, Rn45s abundance, and percent ribosomal gene transcripts. Ribosomal genes were found with the regular expression ^Rp[sl][[:digit:]].
tiss <- ScaleData(object = tiss, vars.to.regress = c("nReads", "percent.ribo","Rn45s"))

Here's a fuller notebook, and we'll have a better organized repository soon.

That's great that the workflow is available. I am just curious why did you decide to use ribosomal genes instead of the mitochondrial genes, which is suggested by Seurat?

Hello! Thanks for asking. In the final version of the paper (this question cites the preprint), we ultimately decided not to normalize on either ribosomal or mitochondrial genes. We found that this normalization added differences where there shouldn't be any, e.g. between cells of the same type in the Bladder.

Thanks for clarifying! I was wondering why I could not find any mention of ribosomal genes anywhere in the methods (in the final paper).


text:
How do PCR duplicates arise and why is it important to remove them for NGS analysis?
I am trying to understand PCR duplicates in NGS analyses (actually whole-genome). I searched, and the best answer I found is in this blog.
However I don't understand if I understood how PCR duplicates arise  correctly because I cannot see the problem of having them in the downstream analysis - aside from computational problems, i.e. unnecessary redundancy.
If I understood correctly, PCR duplicates arise during library preparation when PCR amplifies the fragments with adapters. In this case, if you have duplicate fragments, you'll end up amplifying some of them twice or more.
However, when you are doing the mapping, they should fall on the same region, most likely decreasing the quality of the mapping (since they increase the consensus on a specific sequence, which could have been subject to sequencing errors). But aside from that, you should have the same mapping with respect to the mapping where duplicates are removed, although with decreased quality.
Is the quality problem the actual reason for removing PCR duplicates or is there anything that I'm missing?

In any scenario where depth of coverage is an important factor, PCR duplicates erroneously inflate the coverage and, if not removed, can give the illusion of high confidence when it is not really there.
For example, consider the following hypothetical scenario.
                                         *
TTTCATACTAACTAGCCTGCGGTCTGTGTTTCCCGACTTCTGAGTCATGGGGTTTCAATGCCTATAGATTC
               ..........................C.
                ............................
                      .............T..............
                       ..................C.........
                           ............................
                               ............................
                                   ............................
                                    .....C......................
                                       ..C.........................

The marked position has a coverage of 9 (9 overlapping reads), 4 of which suggest the presence of an alternate allele. Since these 4 reads map to different positions, they are independent observations supporting the same putative single nucleotide variant (SNV). The T in the third read can safely be ignored as a sequencing error, since it only appears once.
Now consider the following hypothetical scenario.
                                         *
TTTCATACTAACTAGCCTGCGGTCTGTGTTTCCCGACTTCTGAGTCATGGGGTTTCAATGCCTATAGATTC
                ............................
                      .............T..............
                       ..................C.........
                       ..................C.........
                       ..................C.........
                       ..................C.........
                               ............................
                               ............................
                               ............................

In this case, the marked position also has a coverage of 9, but 7 of those reads appear to be PCR duplicates. If we remove the duplicates, the marked position has only a coverage of 4, and only one C at the marked position, which is not enough confidence to call the SNV.

Best explanation I've seen so far

I agree with BCArg, very simple and clear explanation. Thank you.

PCR polymerases introduce errors. When an error arises in the first few cycles of amplifications, it will appear in a reasonably high fraction of DNA fragments in the library. After sequencing, you may see the same error occur to multiple reads. If you remove PCR duplicates when calling variants, all errors are reduced down to one read. For high-coverage data, you won't call the error as a false variant. However, if you keep PCR duplicates, the multiple occurrences of the error may make the error look like a real variant. This is particularly bad for calling SNVs from impure tumor samples.


text:
How to wget a gff from NCBI database?
Say I want to use wget to download a gff file for the X chromosome of release 5.57 of the Drosophila melanogaster genome. The NCBI page with the gff file for just this is here.
However, I do not see any way to directly wget the gff from this page. What I can do is click "Send to" at the top, select File, and change the Format dropdown to GFF3, and finally click 'Create File', and the gff file downloads onto my computer.
But then I lose the ability to use the -a flag to log the download details, I can't download the gff using wget from a remote server directly but instead I need to download it locally and then upload it, etc. Is there any way to use wget to download gff files from NCBI programmatically instead of having to download it manually from the website?

Are you looking for these?

@user172818 If you mean my example, not exactly, I'm looking for release a 5.57 gff based on chromosome X, that link lets me download a release 6 whole genome gff. But I think I figured out what to do, and I posted an answer to my question.

my suggestion would be to go via the assembly page.
For your example link, I would

Click the link to the BioProject: PRJNA164
On the right side, under "Related information" choose the assembly
On the assembly page, under "Access the data" on the right side, I would choose "FTP directory for RefSeq assembly"

Once in the RefSeq folder, copy the link to the *_genomic.gff.gz and use wget.
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/215/GCF_000001215.4_Release_6_plus_ISO1_MT/GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.gff.gz

Edit: The initial question requested an older versions, which can be found in RefSeq's archive of prior releases. For this specific Drosophila assembly chromosome X, the version 5 can be found at https://ftp.ncbi.nlm.nih.gov/genomes/archive/old_refseq/Drosophila_melanogaster/RELEASE_5/CHR_X/NC_004354.gff

This does show how to wget from NCBI, although it wouldn't work on the specific case I mentioned. This gives me the whole genome, although I'm curious on getting a single chromosome. Also it will give me release 6 (the latest one), whereas my scenario was for release 5. However I think I found a way to do this.

Pardon, the RefSeq archive seems what you are looking for. Have a look at https://ftp.ncbi.nlm.nih.gov/genomes/archive/old_refseq/Drosophila_melanogaster/RELEASE_5/.  The path to the X chromosome's gff should be straight forward. https://ftp.ncbi.nlm.nih.gov/genomes/archive/old_refseq/Drosophila_melanogaster/RELEASE_5/CHR_X/NC_004354.gff

That's exactly it! Could you edit these links into your main answer so I can mark your comment as accepted?

You're welcome. Glad to have helped

I believe I figured out a way to do this for my particular example (although in general cases Adrian's answer works). To get release 5.57 (an earlier release) of chromosome X of Drosophila melanogaster, I believe you can go to the Datasets Genome page of D. melanogaster, scroll down to the Revision History and select Release 5, and at the top of the page, NCBI lets you click on a 'datasets' or 'curl' box to copy the exact command line argument needed to download a GFF (and other files) of that whole genome. Then, I believe it should be possible to open the gff in Python using a package like gffutils, and filter only for chromosome X rows, and then write those rows to a new gff file.

If you are willing to filter, there is no need for python or writing a whole script. A simple grep -P '^NC_004354.3\t' foo.gff > chrX.gff will do.

Did not realize I could do that, thanks. Although you'd probably have to add a bit to this command to also preserve the header lines. Maybe (grep "^#" foo.gff; grep -P '^NC_004354.3\t' foo.gff) > chrX.gff


text:
Downloading genes from NCBI in fasta format
I'm pretty new in bioinformatics.
I need to download FASTA sequences of several genes. The list of genes I have assembled consists of 140 genes, so I'd rather do this through via code than download each gene manually from the NCBI database.
The genes belong to 2 different organisms (bacteria) I have chosen, most of them are found in both organisms and some are specific to one or the other bacterium.
I've tried running this code which I found here: https://www.ncbi.nlm.nih.gov/books/NBK179288/
esearch -db nucleotide -query "aminoacyl-tRNA\ hydrolase" | efetch -format fasta

The gene of interest being aminoacyl-tRNA hydrolase here. But I don't know how to filter it for the specific organism, Neisseria meningitidis. Tried running this, without luck:
esearch -db nucleotide -query "aminoacyl-tRNA\ hydrolase" | efilter -organism [Neisseria meningitidis] | efetch -format fasta

Appreciate any help :)

Hi how does the list of 140 gene you are interested in look like? do you have gene accessions, or Ids or just gene symbols?, ie symbol: BRCA1 ; id :
672 ; accsession NC_000017.11 , are all the same gene.

what output do you get from your commands? nothing, an error, the wrong gene?

Yes, please [edit] your question and add a few lines of your input file. Presumably, you have a list of genes, so please show us the list so we can understand the input format.

@user12256545 Right now I have the list on a word-document with locus name and gene-ID (that would be "pth" and 61281026 for aminoacyl-tRNA hydrolase). All genes seem to have the same accession nr. so that won't be useful unfortunately.

If I put the species name in quotes, I get expected (by me) output:
% esearch -db nucleotide -query "aminoacyl-tRNA\ hydrolase" | efilter -organism "Neisseria meningitidis" | efetch -format fasta
# outputs:
>NZ_JWNU01000010.1 Neisseria meningitidis strain Nm440716 Scaffold10, whole genome shotgun sequence
GACCGCCTATTCGCTTGCCGCCGGCGGCCCCATCATGCAGGCAGGATTACACGCCTTCACACTCGTACCC
# [snip, tons of output]

Note that filtering down to only N. meningitidis still leaves all of the different sequenced accessions of this species, such that you will still have a lot of different annotated copies of this gene. You will need to filter more specifically to reduce that.

This was helpful, thank you! What would you suggest I use as identifier to filter more specifically? I tried substituting "aminoacyl-tRNA\ hydrolase" with gene-ID (61281026) but that doesn't seem to work.

@Vjon I would suggest trying an NCBI taxid in place of the species name.


text:
How to associate point mutations to quantitative protein levels
My experiment is based on a quantitative continuous variable levels (B12 vitamin) in a sample of patients retrospectively selected based on values that can take (under our criteria):

normal
low levels: under the physiological or normal reference values but not enough to be considered deficiency
deficiency (which is associated with poor nutrition, gastric disease, etc.). These data are excluded.

We are interested in low levels. Our hypothesis is to prove if mutations in 10 genes can be associated with these physiological low levels. These physiological low levels are quantified under the format of maximum and minimum value of B12. The genes can take 3 possible categories:

wt: absence of mutation
het: half of mutation
hom: full mutation

Hypothesis:

I would like, if possible, to try something like a genetic risk score, if multiple het or hom states in multiple genes can be additive. I don't know how to approach this, because there can be multiple permutations, and maybe permutation 1 has a stronger association with low levels than permutation 2.
Any ideas on how to carry this out? I don't really know if mutations are additive or opposite.
Not sure if some clustering approach can be helpful here, like PCA or euclidean distance, linear discriminant analysis... Any advice is welcome.

My idea is to execute a summary of levels of B12 for each gene and see in a boxplot possible differences, perhaps look for some association (or correlation), but I doubt it.
Any idea or previous input or experiment like this would be tremendously welcome.

The solution is supervised machine learning (ML) and is fairly standard analytics for precisely this problem.  PCA or euclidean distance, will not prove this especially distance methods. Discriminant analysis is okay. What you need are complex regression analytics and that broadly includes ML. So IMO  the project needs an analyst. I hope thats helpful and again its just my opinion

Could you recommend me something to begin?

I personally would recommend collaboration, unless someone already published an algorithm specific for purpose. Its a bit like a specialist/expert in military logistics wanting to perform a "commando assault", i.e. it would be really hard, particularly if the key source of help was online.

Yes, you can run an additive model (i.e., B*AA = 2*B*Aa, aa = B*AA - 4*Aa) or a dominance model (B*AA = B*Aa). You can test both of these and see which is a better fit to your data.
However, it is practically known already that the model will be at least predominantly additive.


text:
What command/invocation is used to generate NCBI 16SMicrobial blastdb
I'm looking for the exact invocation used to generate the 16SMicrobial database that you can download from here:
https://ftp.ncbi.nlm.nih.gov/blast/db/
I'm hoping to create the same type of blastdb with the same type of metadata with custom sequences.
Platform isn't an issue, but let's say on ubuntu 14.04 or 16.04.
I would like to replicate the creation of the database as closely as possible.
 The most important feature is the taxonomic information as can be seen here:

The databases on the FTP site contain taxonomic information for each
  sequence, include the identifier indices for lookups, and can be up to
  four times smaller than the FASTA. The original FASTA can be generated
  from the BLAST database using blastdbcmd

Creation of a blastdb using makeblastdb from a set of a fasta sequences is not an issue and can be achieved via:
makeblastdb -in <your_file.fasta> -dbtype nucl -out <database_name>

My question is specifically about the invocation NCBI uses to add the metadata that is present in the NCBI's 16SMicrobial blast database as I am keen to make sure I have replicated the process as closely as possible.

What properties, specifically?

Most crucially the taxonomic information, but I would like to replicate this creation of this database as closely as possible.

Running blastdbcmd -db foo -info provides a little information but I haven't seen anything which will report exactly how a blastdb was created.
A blastdb readme suggests that only the parse_seqids option has been added to the standard parameters.

For those from NCBI, the following makeblastdb commands are recommended:

For nucleotide fasta file:   makeblastdb -in input_db -dbtype nucl -parse_seqids
For protein fasta file:      makeblastdb -in input_db -dbtype prot -parse_seqids

The -parse_seqids requires NCBI fasta headers in the correct format (spec can be found here).

I don't know what specific properties you are considering. As far as I know, that's just a normal blast database, like any other. These are produced using the formatdb command which is part of the ncbi-blast software. You haven't specified what system you are using so I can't help you find and install it, but the command will be something like:
formatdb -p F -i yourSeqs.fasta -n yourDBName

formatdb : the command
-p F : protein: false. This is needed to make databases of nucleotide sequences.
-i yourSeqs.fasta : your input file with all your sequences in fasta format.
-n yourDBName : the name (choose whatever you like) for your database.

You can then use yourDBName as a database to blast against. The above assumes yuou are

I would recommend using the newest ncbi-blast+ suite instead of the (now deprecated?) legacy-blast. In that cases formatdb would become makeblastdb, with slightly different arguments.

@mgalardini that sounds very reasonable. I'm afraid I'm showing my age. This is the only one I've ever used but I last did so more than 10 years ago now. Please post an answer with this new tool you mention. I'm sure it's the Right Way© to go today.

@terdon - whatever works I say ;).  Thanks for the help, I've updated the question to make it (hopefully) more clear.  I am specifically looking to add metadata to a blastdb and looking for a reference or similar for what metadata is added to the 16SMicrobial database.

I am not sure whether you are looking for how to produce taxonomy file or rep_set file (qiime reference file types). But you can produce fasta file (rep_Set) like that:
$blastdbcmd -db 16SMicrobial/16SMicrobial -outfmt '%f' -entry 'all' > 16SMicrobial.fasta

Then you can create QIIME-compatible id_to_tax maps with this code:
https://github.com/mtruglio/QIIME_utilities/blob/master/id_to_tax_mapmaker.sh
Hope this helps.


text:
How is prior.count used by edgeR's cpm
edgeR's cpm function has an argument called prior.count. Based on my understanding of the documentation, it is supposed to be adding a fixed number per sample which is proportional to the library size of said sample. Average of all numbers added to all samples would be equal to prior.count.
However, looking at actual data this does not seem to be the case. Given an imaginary data frame of
df = data.frame(a = c(1,2,3,0),c = c(1,2,3,0)*3, b = c(1,2,3,0)*2)

We can try to calculate log cpms by doing
logCPM = cpm(df,log = TRUE,prior.count = .5)

We can also calculate regular cpms by doing
CPM = cpm(df)

If we wanted to see, what is the number that is added to the CPMs before getting logged, we could do
difference = 2^logCPM - CPM

To make it look pretty lets use tibble
tibble::as.tibble(difference)

# A tibble: 4 x 3
                     a                    c                    b
                 <dbl>                <dbl>                <dbl>
1  25641                25641                25641
2  12821                12821                12821
3 -    0.0000000000582 -    0.0000000000582 -    0.0000000000582
4  38462                38462                38462

Here we see that the number that is added has no relation to the library size of the sample. What I want to learn is, how is this number that is added to each cpm is calculated based on prior.count.
I have been digging through the code but that part goes into C++ territory which hinders my understanding of what is going on

I get different results in the difference table. Which version are you using ? (I am in R 3.5.1 and edgeR 3.22.2)

The prior count ends up getting scaled by the ratio of a library size to the average library size and then multiplied by 2 before getting added to each library size (I'm sure there's a good reason for that, but I don't know what it is). Using the example df data frame in your post, let's walk step-by-step through what cpm() is doing:
df = data.frame(a = c(1,2,3,0),c = c(1,2,3,0)*3, b = c(1,2,3,0)*2)
prior.count = 0.5
# First, we need to calculate a library size
lib.size = colSums(df)
# Calculate the average library size and the adjusted priors
ave.lib = mean(lib.size)
adjusted.prior = prior.count * lib.size / ave.lib
# Update the library sizes
adjusted.lib.size = lib.size + 2*adjusted.prior
# Now we can compute the CPM
customCPM = t((log(t(df) + adjusted.prior) - log(adjusted.lib.size) + log(1000000))/log(2))

The matrix transposition looks really ugly in the last step, but it's needed so R adds things like adjusted.prior to each row of df.
Note that the prior count is only used if you compute log(CPM), it's unused otherwise.

How do you know that edgeR is scaling that way to obtain the adjusted.prior?  By the way sweep  could be used to avoid having to transpose the matrix

It's stated in the source code. I know I can use sweep(), but t() takes fewer characters.

Thanks a lot. By the way was I wrong about the code being within the C++ part? If that R code above was part of the package, I have missed it

You didn't miss it, I translated the C++ code to more accessible R. The C++ code is pretty tough to follow unless you're pretty comfortable with it.

thanks this is an awesome answer and helps me understand prior.count


text:
How can I read FCS files using open source libraries?
FCS is a patented data format used for storing flow cytometry data. The most recent version is FCS3.1. There is some documentation on the format, but there is no information on how to read these files. There are some R packages and a MATLAB code to read an FCS file, but I am looking for standard libraries developed either by the FCS consortium or any other group. I also wish to know if FCS is a subset of an existing standard data format that can be read by a standard library using any programming language.
Finally, I would want to convert these files to an easily readable format like HDF5.

I find strange you ask for standard libraries but you exclude the existing standard tools.

@Llopis the OP isn't rejecting standard tools, on the contrary, he is looking for standard libraries.

@terdon I am not sure how standard libraries are decided but a R/Bioconductor package or a MATLAB one seems to me quite standard (Except if the FCS consortium provides one of course).

@Llopis By standard, I mean is there an OS/language independent interface that allows reading of these files. More importantly, I want to know how FCS files can be read by any general language. I am expecting something like a pseudocode, for this purpose.

I’m still not really sure what you want: the documentation you link to is as close to a general pseudocode as you will get: it contains a detailed, concise specification of the format. That’s the same information that a pseudocode would convey, and it’s the conventional way in which data formats are described (pseudocode is virtually never used for this). That, in addition to the existence of libraries for common analysis languages (R, Python, Matlab), should be sufficient to work with the format.

@WYSIWYG In addition to what Konrad said, the flowCore in  Bioconductor is available for Windows, iOS and Linux.

@KonradRudolph What I want is an instruction on how these files are to be read. The documentation describes the format but I could not figure out the general procedure to read the format which is what I meant by pseudocode. Just assume that I do not want to use bioconductor;  lets say that I want to read these files in C. How should I go about that? Do the FACS machines use a standard mechanism/io library for reading/writing FCS files or they use their own codes? I assumed that if a format is defined by a consortium, they would have provided a standard interface for rw operation.

Why do you say .fcs is a patented format? The disclaimer in the linked documentation just states "Attention is called to the possibility that implementation of this standard may require use of subject matter covered by patent rights".
It may not even be possible to patent it, see: https://patents.stackexchange.com/questions/11968/can-patents-on-computer-file-formats-be-obtained

A few years ago, I wrote a python script to convert FCS files into tab-separated format. It was far from handling all the possibilities that the format description offers, but at least it worked for some of the files produced on one of our machine: http://www.igh.cnrs.fr/equip/Seitz/en_equipe-programmes.html
The format documentation I found enabled decoding (see section 3 of the pdf you mention), but it requires reading data in binary mode.
The general idea of this format (and, I guess, many other binary formats), is that there is a header zone at the beginning of the file with a defined number of fields encoding numbers indicating how the rest of the file is structured. So a first phase is to parse this header, following the description given in the documentation of the format. The information extracted from the header tells where to find the data and how it is encoded, still according to rules described in the format documentation.
In case this may be useful, and for the record, here is the code from the above-mentioned script (after stripping comments, some of which are merely copied from the format documentation, and adding a few ones):
#!/usr/bin/env python
"""This script tries to read FCS flow cytometry data.
Format parsing inspired by information found here:
http://isac-net.org/Resources-for-Cytometrists/Data-Standards/Data-File-Standards/Flow-Cytometry-Data-File-Format-Standards.aspx
"""

import re
# To decode binary-encoded data
import struct
import sys

class Parameter(object):
    """This object represents one of the parameter types that are present in a DATA segment of a FCS file."""
    __slots__ = ("p_name", "p_bits", "p_range", "p_ampl", "parser")
    def __init__(self, p_name, p_bits, p_range, p_ampl):
        self.p_name = p_name
        self.p_bits = p_bits
        self.p_range = p_range
        self.p_ampl = p_ampl
        # Function for parsing a value of the parameter in the data segment
        self.parser = None

##############################################
# Here starts the parsing of the header part #
# which tells where the other parts are.     #
##############################################

f = open(sys.argv[1], "rb")
# The format name is encoded in 6 letters
# An ASCII letter is coded with one octet
file_format = "".join([f.read(1) for __ in range(6)])
sys.stdout.write("Format: %s\n" % file_format)
# The format descriptions reserves 4 octets that we skip
skip = f.read(4)
# 8 octet chunks encode the start and end positions
# of different parts of the data
text_start = int(f.read(8).strip(" "))
text_end = int(f.read(8).strip(" "))
data_start = int(f.read(8).strip(" "))
data_end = int(f.read(8).strip(" "))
analysis_start = int(f.read(8).strip(" "))
analysis_end = int(f.read(8).strip(" "))

if (analysis_start and analysis_end):
    sys.stderr.write("Cannot deal with ANALYSIS segment of an FCS file.\n")

####################################################
# Here starts the parsing of the "TEXT" portion    #
# which describes how the data proper is organized #
####################################################
f.seek(text_start)
# The first character in the primary TEXT segment is the ASCII delimiter character.
sep = f.read(1)
if sep not in ["_", "@"]:
    alt_sep = "_@_"
elif sep not in ["_", "|"]:
    alt_sep = "_|_"
else:
    assert sep not in ["+", "|"]
    alt_sep = "+|+"
text_segment = f.read(text_end - text_start)

fields = text_segment.split(sep)

info = {}

i = 0
while i < len(fields) - 1:
    key = fields[i]
    i += 1
    val = fields[i]
    i += 1
    # Keywords are case insensitive, they may be written in a file in lower case, upper case, or a
    # mixture of the two. However, an FCS file reader must ignore keyword case. A keyword value may
    # be in lower case, upper case or a mixture of the two. Keyword values are case sensitive.
    info[key.upper()] = val
print "%s events were detected." % info["$TOT"]
print "Each event is characterized by %s parameters" % info["$PAR"]

if info["$NEXTDATA"] != "0":
    sys.stderr.write("Some other data exist in the file but hasn't been parsed.\n")

# L - List mode. For each event, the value of each parameter is stored in the order in which the
# parameters are described. The number of bits reserved for parameter 1 is described using the
# $P1B keyword. There can be only one set of list mode data per data set. The $DATATYPE
# keyword describes the data format. This is the most versatile mode for the storage of flow
# cytometry data because mode C and mode U data can be created from mode L data.
assert info["$MODE"] == "L"

parameters = []

# indices of the parameters
p_indices = range(1, int(info["$PAR"]) + 1)
for i in p_indices:
    p_name = info["$P%dN" % i]
    p_bits =  info["$P%dB" % i]
    p_range = info["$P%dR" % i]
    p_ampl =  info["$P%dE" % i]
    parameters.append(Parameter(p_name, p_bits, p_range, p_ampl))

sys.stdout.write("The parameters are:\n%s\n" % "\t".join([par.p_name for par in parameters]))

# How are 32 bit words organized
if info["$BYTEORD"] == "4,3,2,1":
    endianness = ">"
else:
    endianness = "<"
    assert info["$BYTEORD"] == "1,2,3,4"

    # I stripped a long comment which is just a copy of the documentation
# Type of data:
if info["$DATATYPE"] == "I":
    for par in parameters:
        nb_bits = int(par.p_bits)
        assert nb_bits % 8 == 0
        nb_bytes = nb_bits / 8
        # Determine format string for unpacking (see https://docs.python.org/2/library/struct.html)
        if nb_bytes == 1:
            c_type = "B" # unsigned char
        elif nb_bytes == 2:
            c_type = "H" # unsigned short
        elif nb_bytes == 4:
            c_type = "L" # unsigned long
        elif nb_bytes == 8:
            c_type = "Q" # unsigned long long
        else:
            raise ValueError, "Number of bytes (%d) not valid for an integer (see https://docs.python.org/2/library/struct.html#byte-order-size-and-alignment)." % nb_bytes
        fmt = "%s%s" % (endianness, c_type)
        p_range = int(par.p_range)
        def parser(data):
            value = struct.unpack(fmt, data.read(nb_bytes))[0]
            try:
                assert value < p_range
            except AssertionError:
                print "Value %s higher than %d" % (str(value), p_range)
            return value
        par.parser = parser
    pass
else:
    raise NotImplementedError, "Only the parsing of integer value has been implemented so far."

out_file = open(sys.argv[2], "w")
out_file.write("#amplification_types\t" + "\t".join([par.p_ampl for par in parameters]) + "\n")
out_file.write("parameters\t" + "\t".join([par.p_name for par in parameters]) + "\n")
i = 1
##############################################
# Here starts the parsing of the data proper #
##############################################
f.seek(data_start)
while f.tell() < data_end:
    values = []
    for par in parameters:
        values.append(par.parser(f))
    out_file.write("%d\t" % i + "\t".join(map(str, values)) + "\n")
    i += 1
out_file.close()
f.close()

The general idea of this format (and, I guess, many other binary formats), is that there is a header zone at the beginning of the file with a defined number of fields encoding numbers indicating how the rest of the file is structured. So I begin by parsing this header, following the description given in the documentation of the format. The information extracted from the header tells me where to find the data and how it is encoded, still according to rules described in the format documentation. I'm not good at writing pseudocode, but I can try to add comments in the python code.

I maintain several Python libraries for interacting with flow cytometry data (FCS files). The most recent of which, FlowKit, provides functionality to export the event data as a NumPy array, Pandas DataFrame, or CSV text. Using the h5py library, it's trivial to export as an HDF5 file

R/Bioconductor has a number of different flow cytometry processing packages. One place to start for looking at cytometry data from a high level would be openCyto (or its vignette), which is a large set of tools for basic extraction and analysis of FCS files.
I have looked in the past at the FCS files as an R structure using flowCore. Loading a single FCS file is fairly straightforward and follows a familiar R pattern:
file.name <- "/dir/file.fcs"
x <- read.FCS(file.name, transformation=FALSE)
summary(x)

Asking about converting FCS files into an "easily readable format like HDF5" doesn't seem like the right question. HDF5 is a container format and shares a lot of similarity to file systems. I've found it best to keep FCS files as they are, as it is a compact, standardised binary format.

HDF5 is not a major concern (I just like it). I am still curious about how something like read.FCS works. I believe if FCS is a standard format then there should be standard (non-proprietary) libraries to read/write FCS files. As of now I'm forced to use R or Matlab (and still there are no standard libraries here too. I guess there are just different packages/programs made by different groups/individuals; no single standard). What if I want to analyse a FCS file on C, Fortran or Python (etc)?


text:
How to retrieve the yeast locus tag from uniprot via sparql?
Let's say I want to retrieve the yeast-specific identifier for a certain protein, in the example below it will be P00330 which I would like to link to YOL086C.
When I go to uniprot's sparql UI and type
PREFIX up:<http://purl.uniprot.org/core/>
PREFIX rdfs:<http://www.w3.org/2000/01/rdf-schema#>

SELECT *
WHERE
{
    Bind(<http://purl.uniprot.org/uniprot/P00330> as ?protein).
    ?protein rdfs:seeAlso ?db .
    ?db up:database <http://purl.uniprot.org/database/SGD> .

}

it returns
protein db
http://purl.uniprot.org/uniprot/P00330  http://purl.uniprot.org/sgd/S000005446

How would I now go on and find the YOL086C? When I add the line
?db ?p ?o

to the query above, I don't see anything that would directly lead me to the identifier.
I also tried to access another database
?db2 up:database <http://purl.uniprot.org/database/EuPathDB>

but this returns all kind of databases while I am interested only in the FungiDB. How could I select this one?
How could I link the two IDs to each other?

Ok, figured it out by looking at the scheme uniprot provides:
PREFIX up:<http://purl.uniprot.org/core/>

SELECT ?protein ?yeast_id
WHERE
{
    Bind(<http://purl.uniprot.org/uniprot/P00330> as ?protein).
    ?protein up:encodedBy ?gene .
    ?gene up:locusName ?yeast_id
}

So one starts at the accession number in the middle and then uses the path to the actual locusName.
Also the path via the database works fine, I just did not use the seeAlso correctly (protein was only bound to db but not to db2, that's why it returned all entries):
PREFIX up:<http://purl.uniprot.org/core/>
PREFIX rdfs:<http://www.w3.org/2000/01/rdf-schema#>

SELECT ?protein ?db ?db2
WHERE
{
    Bind(<http://purl.uniprot.org/uniprot/P00330> as ?protein).
    ?protein rdfs:seeAlso ?db .
    ?protein rdfs:seeAlso ?db2 .
    ?db up:database <http://purl.uniprot.org/database/SGD> .
    ?db2 up:database <http://purl.uniprot.org/database/EuPathDB>

}

This will return the desired output (one would still have to process the db2 entry, of course):
protein db  db2
http://purl.uniprot.org/uniprot/P00330  http://purl.uniprot.org/sgd/S000005446  http://purl.uniprot.org/eupathdb/FungiDB:YOL086C

I am really happy to see you found your answer! Could we use this as an example query on the sparql endpoint itself?

@Jerven: feel free to use it, sure! I just updated my answer to show that also the path via the database works fine; I just used it wrong.


text:
How to retrieve InChI key for KEGG compound?
I would like to retrieve the InChI representation for a given KEGG compound but I fail to find a direct solution for this (preferably doing it via bioservices).
One could do it via ChEBI like this:
from bioservices import *

kegg_con = KEGG()

kegg_entry = kegg_con.parse(kegg_con.get('C00047'))

chebi_con = ChEBI()

chebi_entry = chebi_con.getCompleteEntity('CHEBI:' + kegg_entry['DBLINKS']['ChEBI'])

print chebi_entry.smiles
print chebi_entry.inchi
print chebi_entry.inchiKey

which will print
NCCCC[C@H](N)C(O)=O
InChI=1S/C6H14N2O2/c7-4-2-1-3-5(8)6(9)10/h5H,1-4,7-8H2,(H,9,10)/t5-/m0/s1
KDXKERNSBIXSRK-YFKPBYRVSA-N

Is there a direct way which does not require to go via ChEBI?

I had a look at the KEGG rest API but couldn't find anything; I think you might need to access the data in the KEGG FTP ($$$)

@mgalardini: I checked that, but seems it is not available at all, or I overlooked it. The way above is still feasible, just requires one additional step.

It's not too bad, especially considering that you have a library to do it for you :)

I'm not aware of a functionality in the KEGG service that would allow you to do that directly. I believe your solution (using e.g. ChEBI) is the correct one.
You may have other useful service such as unichem, chemspider (you would need a login though for that one) could help as well.
Disclaimer: I'm the bioservices main author


text:
How to access organism IDS, taxon nodes from UniProt?
I try to find organism IDs used by UniProt (e.g. 9606 for "Homo sapiens") and information about the 'taxonomic kingdom' this organism belongs to. The only source on UniProt I could find is this file which contains all required information. So in the Python world it would be a dictionary of the form:
{org_id: {'name': 'org_name', 'kingdom': 'org_kingdom'}}

or a dataframe.
However, the file from above is not particular nice to parse due to its inconsistent column format (not impossible to parse though), so I am wondering whether there is an alternative source to retrieve this information or whether there is maybe even a webservice available that takes care of this parsing.

That file is actually trivial to parse. You only care about lines that have N=$species, so you can simply do:
sed -En 's/.* ([0-9]+): N=(.*)/\1\t\2/p' speclist

That will return a tab separated list of taxID and species:
$ sed -En 's/.* ([0-9]+): N=(.*)/\1\t\2/p' speclist | head
648330  Aedes albopictus densovirus (isolate Boublik/1994)
10804   Adeno-associated virus 2
648242  Adeno-associated virus 2 (isolate Srivastava/1982)
118452  Abacion magnum
72259   Abaeis nicippe
102642  Abax parallelepipedus
392897  Abalistes stellaris
75332   Abbottina rivularis
515833  Abdopus aculeatus
56673   Antarctic bacterium DS2-3R

So to get human, for example, you could do:
$ sed -En 's/.* ([0-9]+): N=(.*)/\1\t\2/p' speclist | grep sapiens
63221   Homo sapiens neanderthalensis
9606    Homo sapiens

And, for exact matches:
$ sed -En 's/.* ([0-9]+): N=(.*)/\1\t\2/p' speclist |
    awk -F'\t' '$2=="Homo sapiens"'
9606    Homo sapiens

If you also want to keep the taxon node (kingdom), you can use:
sed -En 's/.* ([A-Z])  *([0-9]+): N=(.*)/\1\t\2\t\3/p' speclist

Which prints:
$ sed -En 's/.* ([A-Z])  *([0-9]+): N=(.*)/\1\t\2\t\3/p' speclist  | head
V   648330  Aedes albopictus densovirus (isolate Boublik/1994)
V   10804   Adeno-associated virus 2
V   648242  Adeno-associated virus 2 (isolate Srivastava/1982)
E   118452  Abacion magnum
E   72259   Abaeis nicippe
E   102642  Abax parallelepipedus
E   392897  Abalistes stellaris
E   75332   Abbottina rivularis
E   515833  Abdopus aculeatus
B   56673   Antarctic bacterium DS2-3R

The sed command might appear a little daunting if you don't know sed, but it's actually quite simple. The -E flag enables extended regular expressions which give us a simpler syntax. The -n suppresses normal output, so that sed won't actually print anything unless explicitly told to.
The s/from/to/ is the substitution operator and will replace from with to. Here, we are replacing zero or more characters (.*) until a space followed by one or more numbers ([0-9]+), a :, another space, then the string N= and then everything else (.*). The parentheses let us capture the matched strings, so that we can then refer to the 1st captured string as \1, the second as \2 etc. We are therefore replacing the entire line (because the regex matches the entire line) with the 1st pattern captured (the taxID, \1), a tab (\t) and the species name (\2). Finally, the p at the end means "print the resulting line if the substitution was successful", so only lines matching the desired pattern are actually printed.

Thanks, that helps already. Is there a quick way to also extract the taxonomic kingdom and store it in the third column?

@Cleb sure, see updated answer. Next time, please explain all these details in your question (I mean tell us exactly what information you want to extract) so we can give you a working answer directly.

These are NCBI taxonomy ID's rather than anything specific to UniProt. One way to parse them would be using BioPython via Bio.Entrez, for example:
from Bio import Entrez

def get_taxon_name_from_NCBI_id(id_):
    Entrez.email = 'name@provider.com'
    handle = Entrez.efetch(db='taxonomy', id= id_, retmode='xml')
    records = Entrez.read(handle)[0]
    return records['ScientificName'], records['TaxId']

print(get_taxon_name_from_NCBI_id(9606))
# ('Homo sapiens', '9606')

Of course you could query Entrez directly:
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=taxonomy&id=9606
Or even more basically do a manual search of NCBI taxonomy:
https://www.ncbi.nlm.nih.gov/taxonomy/?term=9606

I was not precise enough in my question: What I am after are the IDs, so basically what you use as input and there associated name and taxonomic kingdom. Ideally it would be a dictionary of the form {org_id: {'name': 'org_name', 'kingdom': 'org_kingdom'}} but I guess I can easily do this directly from NCBI; thanks for making me aware of this!


text:
How to extract the protein fasta file from a genbank file?
Let's say I have a genbank file, e.g. the one downloaded from here, that contains entries as
 gene            190..255
                 /gene="thrL"
                 /locus_tag="b0001"
                 /gene_synonym="ECK0001; JW4367"
                 /db_xref="EcoGene:EG11277"
                 /db_xref="GeneID:944742"
 CDS             190..255
                 /gene="thrL"
                 /locus_tag="b0001"
                 /gene_synonym="ECK0001; JW4367"
                 /function="leader; Amino acid biosynthesis: Threonine"
                 /GO_process="GO:0009088 - threonine biosynthetic process"
                 /codon_start=1
                 /transl_table=11
                 /product="thr operon leader peptide"
                 /protein_id="NP_414542.1"
                 /db_xref="ASAP:ABE-0000006"
                 /db_xref="UniProtKB/Swiss-Prot:P0AD86"
                 /db_xref="EcoGene:EG11277"
                 /db_xref="GeneID:944742"
                 /translation="MKRISTTITTTITITTGNGAG"

and I want to extract the potein fasta from it, how would I do this?
I found a code sample here. When I then run:
from Bio import SeqIO

gbk_filename = 'NC_000913.3.gb'
faa_filename = 'NC_000913.3_warwick.fasta'
input_handle  = open(gbk_filename, "r")
output_handle = open(faa_filename, "w")

for seq_record in SeqIO.parse(input_handle, "genbank") :
    print "Dealing with GenBank record %s" % seq_record.id
    for seq_feature in seq_record.features :
        if seq_feature.type=="CDS" :
            assert len(seq_feature.qualifiers['translation'])==1
            output_handle.write(">%s from %s\n%s\n" % (
                   seq_feature.qualifiers['locus_tag'][0],
                   seq_record.name,
                   seq_feature.qualifiers['translation'][0]))

output_handle.close()
input_handle.close()

it results in
     31     for seq_feature in seq_record.features :
     32         if seq_feature.type=="CDS" :
---> 33             assert len(seq_feature.qualifiers['translation'])==1
     34             output_handle.write(">%s from %s\n%s\n" % (
     35                    seq_feature.qualifiers['locus_tag'][0],

KeyError: 'translation'

How could this be fixed?

One can get it to work by using SeqIO.InsdcIO.GenBankCdsFeatureIterator:
from Bio import SeqIO

file_name = 'NC_000913.3.gb'

# stores all the CDS entries
all_entries = []

with open(file_name, 'r') as GBFile:

    GBcds = SeqIO.InsdcIO.GenBankCdsFeatureIterator(GBFile)

    for cds in GBcds:
        if cds.seq is not None:
            cds.id = cds.name
            cds.description = ''
            all_entries.append(cds)

# write file
SeqIO.write(all_entries, '{}.fasta'.format(file_name[:-3]), 'fasta')

This will then give the desired output (just the first two entries):
>b0001
MKRISTTITTTITITTGNGAG
>b0002
MRVLKFGGTSVANAERFLRVADILESNARQGQVATVLSAPAKITNHLVAMIEKTISGQDA
LPNISDAERIFAELLTGLAAAQPGFPLAQLKTFVDQEFAQIKHVLHGISLLGQCPDSINA
ALICRGEKMSIAIMAGVLEARGHNVTVIDPVEKLLAVGHYLESTVDIAESTRRIAASRIP
ADHMVLMAGFTAGNEKGELVVLGRNGSDYSAAVLAACLRADCCEIWTDVDGVYTCDPRQV
PDARLLKSMSYQEAMELSYFGAKVLHPRTITPIAQFQIPCLIKNTGNPQAPGTLIGASRD
EDELPVKGISNLNNMAMFSVSGPGMKGMVGMAARVFAAMSRARISVVLITQSSSEYSISF
CVPQSDCVRAERAMQEEFYLELKEGLLEPLAVTERLAIISVVGDGMRTLRGISAKFFAAL
ARANINIVAIAQGSSERSISVVVNNDDATTGVRVTHQMLFNTDQVIEVFVIGVGGVGGAL
LEQLKRQQSWLKNKHIDLRVCGVANSKALLTNVHGLNLENWQEELAQAKEPFNLGRLIRL
VKEYHLLNPVIVDCTSSQAVADQYADFLREGFHVVTPNKKANTSSMDYYHQLRYAAEKSR
RKFLYDTNVGAGLPVIENLQNLLNAGDELMKFSGILSGSLSYIFGKLDEGMSFSEATTLA
REMGYTEPDPRDDLSGMDVARKLLILARETGRELELADIEIEPVLPAEFNAEGDVAAFMA
NLSQLDDLFAARVAKARDEGKVLRYVGNIDEDGVCRVKIAEVDGNDPLFKVKNGENALAF
YSHYYQPLPLVLRGYGAGNDVTAAGVFADLLRTLSWKLGV

Is there any way to include the proteins names annotated at the gbk file as the sequence description?