text:
Is there a tool to convert a variableStep .wig to fixedStep one?
I have a .wig file in this format:
variableStep chrom=chr1 span=20
59965   7
59985   10
60005   10
60025   10

and I would like to format it to something like this format:
#bedGraph section chr16:0-103584
chr16   0       59901   0
chr16   59901   59927   1
chr16   59927   59931   2

Please assume, I am working with a reference hg19 genome, in case the info isn't deducible from the variableStep .wig file. I couldn't find anything in WiggleTools

To me it is not obvious if these two files are consistent. What is the first position 59965 converted to ? (Perhaps it is because I'm unfamiliar with these type of files..)

@llrs they aren’t of the same data. It’s just to demonstrate the format

Could you please change one or the other so that they represent the same data?

When I had to do something similar previously, I had used the following method using sed and awk. It is not as robust in terms of error detection and correction as the python solution but it is a simple one-liner.
$ cat temp.wig
variableStep chrom=chr1 span=20
59965   7
59985   10
60005   10
60025   10

$ sed -r 's/variableStep chrom=(.*) span=([0-9]*)/\1\t\2/g' temp.wig \
        | awk 'BEGIN{FS="\t";OFS="\t"}{if ($1~/^chr/) {c=$1 ; s=$2} else {print c,$1,$1+s,$2}}' \
    > temp.bedgraph

$ cat temp.bedgraph
chr1    59965   59985   7
chr1    59985   60005   10
chr1    60005   60025   10
chr1    60025   60045   10

This is a temporary solution I came up with that worked for my specific case:
def convert_line(line, chrom, span, start=None, fixedStep=None, count=None):
        new_line=''
        if fixedStep is None or count is None or start is None:
                line_arr  = line.split()
                start = line_arr[0]
                end   = int(line_arr[0]) + int(span)
                val   = line_arr[1]
                new_line = '{0}\t{1}\t{2}\t{3}\n'.format(chrom, start, end, val)
        else:
                line_arr  = line.split()
                val   = np.float(line_arr[0])
                start = int(start) + int(fixedStep)*int(count)
                end   = start + int(fixedStep)
                new_line = '{0}\t{1}\t{2}\t{3}\n'.format(chrom, start, end, val)

        return new_line

def convert_file(input_file, output_file):
        chrom=None
        span=None
        fixedStep=None
        count=None
        start=None
        with open(input_file, 'r') as fin, open(output_file, 'w') as fout:
                for line in fin:
                        if len(line)<1 or line[0] == '#':
                                continue
                        elif 'variableStep' in line or 'fixedStep' in line:
                                for key in line.split():
                                        if 'chrom=' in key:
                                                chrom = key.split('chrom=')[-1]
                                        if 'span=' in key:
                                                span = key.split('span=')[-1]
                                        if 'step=' in key:
                                                fixedStep = key.split('step=')[-1]
                                        if 'start=' in key:
                                                start = key.split('start=')[-1]
                                                count = 0
                        else:
                                line_out = convert_line(line, chrom, span, start, fixedStep, count)
                                fout.write(line_out)
                                if count is not None:
                                        count += 1

input:
variableStep chrom=chr16 span=20
59965   7
59985   10
60005   10

output:
chr16   59965   59985   7
chr16   59985   60005   10
chr16   60005   60025   10
chr16   60025   60045   10


text:
How to use analysis output in my pipeline
I am running an analysis and I want to edit the output that I received from that analysis. However I cannot define the output name in my pipeline, because after running the analysis it is asking me to name the output file. How can I include that output file in my pipeline ?

Welcome to Bioinformatics.StackExchange! I would recommend reverting the edit you just made. The question is impossible to answer generically, since many tools have different ways to specify output filenames. My answer is only applicable to pepstats, like your original question asked.

I assume you are referring to the pepstats program from EMBOSS. The program has a -outfile flag that allows you to specify your desired name for the output file when you're not running the program interactively.

source: Stack Exchange

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

sort: votes


text:
bedtools: get name column of alignment file A
I'm using the coverage function of bedtools to check whether a certain set of regions (file B) has any overlap with known regions of epigenetic markers (file A). When that is the case, I want to access the name of that marker, which is given as the third column (name) in file A.
However, when I run bedtools coverage -a fileA.bed -b fileB.bed, I only get the columns of file B as an output, together with the overlap fraction. Is there a built-in way of accessing also the columns of the file A to which I aligned, or an easy way of cross-referencing to it?

I would recommend using bedtools intersect for that type of comparison.  This will allow you to see the specific segments that overlap.  bedtools coverage is used more for determining how many things cover a region but not which ones.
The documentation linked above will provide a much better desciption of intersect than I would be able to in this post but I would recommend using the -loj or -u depending on what you want.
-loj    Perform a “left outer join”. That is, for each feature in A report each overlap with B. If no overlaps are found, report a NULL feature for B.

-u  Write original A entry once if any overlaps found in B. In other words, just report the fact at least one overlap was found in B. Restricted by -f and -r.

If you're working with fileB.bed (a sorted BED file) and fileA.txt (a three-column text file, where the name is in the third column), then you could pipe single-base pseudo-BED elements to bedmap.
$ bedmap --echo --echo-map-id --delim '\t' fileB.bed <(awk -vFS="\t" -vOFS="\t" '{ print $1, $2, ($2+1), $3 }' fileA.txt) > answer.bed

The --echo-map-id option reports the list of marker IDs from fileA.txt, where the single-base marker interval from fileA.txt overlaps a region in fileB.bed.
If all the marker IDs are not unique (if there are duplicate marker IDs) it may make sense to instead use --echo-map-id-uniq.
This assumes that marker positions in the second column and that the intervals in fileB.bed are using a 0-based index scheme. If marker positions are 1-based and fileB.bed is 0-based, then a small adjustment is needed to make both inputs use the same indexing, to calculate correct overlaps:
$ bedmap --echo --echo-map-id --delim '\t' fileB.bed <(awk -vFS="\t" -vOFS="\t" '{ print $1, ($2-1), $2, $3 }' fileA.txt) > answer.bed

I mention this, as SNVs (such as those from VCF) can often be reported with a 1-based index scheme, which adds some nuance to doing set operations against BED files.

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


text:
I need some tips and suggestions for further analysis of NGS expression data (log2cpm)
I am a PhD student who inherited some log2cpm data of expression data from bulk kidney tissue from a UUO(unilateral urethral obstruction) experiment that tests a new drug. The sample material consists of:

6 x Ligated Kidney (Untreated
6 x Ligated Kideny (Treated)
3 x Unligated Kidney (Untreated)
3 x Unligated Kidney (Treated)

The previous study showed that this drug was effective against fribrosis, so the aim of my study is to investigate how this drug affects inflammation and mitochondrial function in the diseased kidney.
So I have filtered my original dataset based on MitoCarta v2 DB, so investigate how the drug affected mitochondrial function. I afterwards run a PCA analysis on log2cpm data which revealed the following:

The following figure shows PC1(84.51%) and PC2(4.69%).
Red (Ligated, Untreated), purple (Ligated, Treated), green(healthy kidney, not treated), orange (Healthy kidney, treated).
I interpret this as there is a very small difference between treated and not treated ligated kidneys, however there is a small difference.
I have just started doing bioinformatics and I am very unsure which approach I should choose further, however I have though maybe do a highly variable genes analysis (HVG) see link, isolate the 500 most variable genes and run a heatmap with hieraichal clustering. I was also thinking of looking at the loadings in PC2, and take this further. The truth is that I am very unsure what is the right approach, and my research group does not have bioinformatician I can ask anymore.
I have the following questions:

What kind of analysis you guys would suggest I use to investigate
further (Hieraichal clustering, HVG etc?)
Is there any alternative preprocessing I could use before I run the PCA, to make the analysis more accurate?
Are there any packages that you can suggest?
Any other tips?

Any other tips?
I asked this same question on Biostars as well
https://www.biostars.org/p/367722/#367730

Welcome to the site. What is your expected result? A list of genes to then validate experimentally with an immuno* or some list of genes for prediction of status/response... What other information do you know about those samples: Could it be that some samples are before and after treatment? Have you considered analysing the differentially expressed genes on the whole dataset (without filtering first for mitochondrial related genes)?

It looks like you'd be best served with a standard differential expression analysis using a factorial design. You should get the raw counts data for that though.

Hey,
the  Animals were sacrificed on day 7 or day 15 after the surgery, and treatment started immediatley after surgery. Drug was given immediatley after surgery.
The reason I want to investigate mitochondrial and immune genes, is that I previously did a GSEA with TopGO package of the whole transcriptome based on logfc-values and q-values (down and upregulated genes isolated, qval-cutoff 0.01, (+/-)fc 1.5) of Ligated(untreated vs. healthy kidneys), and Ligated(treated) vs. Ligated(untreated) kidneys. see figure. It showed that mosltey mitochondrial and inflammatin was affected

Now I want to investigate a bit further what molecular pathways the drug affects and what constitute the improvement.

@NewbieCoder Did you include the time in the design matrix when analysing differential expressed genes? Also if you are interested in those pathways is different than being interested in the genes expressed in the mithocondria as it appeared first. Could you [edit] the question and clarify what is your question?

Hey, I can add some information tomorrow. I was not at work the past days :) Then I can clarify exactly what I need to do. If you have any tips, I would appreciate it greatly :)

Since you indicate already having performed differential expression and using gene set enrichment the next step is to look at what changes due to drug treatment in the pathways of interest while having a decent understanding of those pathways. Things like IPA or pathview (from Bioconductor) can aid in this, but at the end of the day you need to have a decent understanding of the pathways yourself to judge what seems to be like the best targets for future work.


text:
Error while loading a big file with Cytoscape
I installed Cytoscape 3.7.1 with java 1.8.0_191 in windows server 2012. I have a 7.4 GB csv file (about 1,500,000 reords) and when I tried to load it into the Cytoscape it throws an error java.lang.OutOfMemoryErroe:Java heap space.
My system has 500GB ram and then I changed Cytoscape.vmoptions file according to -Xms102400M but when I load my file, in task manager, the memory increased until 50GB and then Cytoscape crashed and again it threw the same error: java.lang.OutOfMemoryErroe:Java heap space.
My data is like this sample:
45677888,test1,3453453453,test2,3534523453,235412352,3453452345,test3,235423452,test4,1980/09/02,14:13,23523525,test5
45234288,test11,34234553453,test12,353434534553,23453452352,342422345,test13,23456543452,test14,1980/04/01,14:12,2323423425,test15
4243424888,test12,3235253453,test22,3533456343453,233534532352,2343452345,test33,23345342,test44,1980/03/4,14:11,23674575,test55

My file is in CSV format. col2 is source node, col1 and col3 are attributes for source node,col 4 is destination node, col5,col6,col7 are attributes for destination node, col8 until end are attributes for edge. I have about 200,000 nodes and 230000 edges.
How can I read this file in Cytoscape?

How much RAM does your computer have? And how big is this network? 7.4G csv? Can you show us a few lines? Does it have a lot of extra data or is it just a list of edges?

thank you for your response. My answers to your questions add to my post.

Please don't post images of text! Just copy the text and paste it directly into your question. Use the formatting tools to format it as code. And could you explain what format this is? Where's the node? Where's the edge? How many edges and nodes do you have? And, more importantly, why do you want to load such an enormous network into a vizualization tool like cytoscape? You will probably not be able to do anything useful with such a large file anyway. Can't you limit the size somehow?

this data is a connected network and its important for me that I can see it as much as possible on one page.

You won't be able to see anything useful in a GUI with a network of this size. Even if you do manage to load it.

if i first detect community in my network and assign a node for each community and then by click on each community, i can see all nodes of that community. for this, what solution do you recomment?

On a network of 200,000 nodes? Even if cytoscape manages to somehow deal with something that huge, I really doubt you will be able to visualize it in any useful way. f you're looking for communities, why don't you first extract the communities into subgraphs and load those into cytoscape? Alternatively, can you isolate only the connected component of the network? That should also reduce the size. 23e4 edges for 20e4 nodes is a very sparse network. You should be able to separate that into smaller, connected sub networks. There's no point in looking at the whole thing if it isn't connected.

thanks for your time. My goal is to review a connected network with 200,000 nodes in cytoscape. To reach this goal and for the greate of the network, I decide first detect community and show each community as a node that edges between this nodes is edges between communities. Then by click on each node, I can see relation between nodes in a community. Can I do this solution by cytoscape?

I really don't think you will be able to do anything useful in cytoscape with such a huge network. As you can see, it can't even read it. Imagine how hard it will be for cytoscape to manipulate it! Why don't you just use the communities as nodes? Is that an option for you?

what tool do you recommend? i can not just use the communities as nodes. relation between nodes in communities are important for me.

I don't think any graphical tool will be able to handle this amount of data. I would instead try two things: i) make a different network for each community and use that to analyze the relationships of the nodes within the community ii) make another network where the communities are the nodes and use that to analyze the relationship between the communities.

The first thing you can try is to separate the network from the meta information. Assuming you have access to a *nix machine, you can run:
awk -F"," '{print $2,$4}' file.csv > network.csv

That will produce a file like this (based on your example data):
test1 test2
test11 test12
test12 test22

That should radically decrease the file size and will give you the best chance of loading it in Cytoscape. You can then similarly extract the attributes:
awk -F"," '{print $2,$1}' file.csv > source.column1.attributes
awk -F"," '{print $2,$3}' file.csv > source.column3.attributes
awk -F"," '{print $4,$5}' file.csv > target.column5.attributes
awk -F"," '{print $4,$6}' file.csv > target.column6.attributes
awk -F"," '{print $4,$7}' file.csv > target.column7.attributes
awk -F',' '{print $2,$4,$8}' file.csv > edge.column8.attributes
awk -F',' '{print $2,$4,$9}' file.csv > edge.column9.attributes
awk -F',' '{print $2,$4,$10}' file.csv > edge.column10.attributes
awk -F',' '{print $2,$4,$11}' file.csv > edge.column11.attributes
awk -F',' '{print $2,$4,$12}' file.csv > edge.column12.attributes
awk -F',' '{print $2,$4,$13}' file.csv > edge.column13.attributes
awk -F',' '{print $2,$4,$14}' file.csv > edge.column14.attributes

Now, try loading network.csv into Cytoscape and then load each of the node and edge attributes separately.

However, a network with 200,000 nodes and 230,000 edges is just too large to be able to manipulate in any useful way in a GUI. It will just look like a horrible hairball. I really urge you to do some preprocessing and pare the network down to something more manageable. I used to work with a human PPI network of ~13000 nodes and ~15000 edges. While I could load that into Cytoscape, that was already far too big to be able to do anything useful with. I could use the data analysis aspects of Cytoscape, but not really the visualization.

I've done some testing and ~50,000 nodes with ~50,000 edges is the upper limit of what my 2015 macbook pro can handle. Even then it is painfully slow though. Anything you can do to reduce to graph will help.

@conchoecia well, the OP is using a machine with 500GB, presumably more than your laptop :). So I am guessing their limit will be considerably higher. But even so, even if they do manage to load this behemoth of a network onto Cytoscape, you can't really do anything useful with a GUI on something so large.

thanks all, that is right. if i can load such enormous network into cytoscape, GUI of it is not useful, but unfortunately, cytoscape displays data as soon as it is loaded.it is better that it load data without display and then we can filter on data and then display network.

I work a lot with large networks and I agree with what most other posters have said, in that it is hardly useful to explore a network like this visually and Cytoscape - while it is one of best GUIs in general - is pretty bad with large graphs. Here are some alternatives you could try if you want to go down that road:

Gephi also probably won't handle that number of nodes/edges well, but it may not crash. Development on it has stagnated though.
Graphia I haven't tried this, but it is supposed to be better at handling large graphs. This project is in its infancy, so don't expect too much

If you care about the analysis of the graph itself, I'd go with NetworKit or SNAP which are much better suited to large graphs.
If you want to visualize the networks, something that tends to perform well here is:

Graph embeddings with VERSE
Dimensionality reduction of the embeddings with UMAP
Plot with Datashader

This can handle absolutely massive graphs if you have the compute - which you seem to have. Here's one I am just working on (albeit small compared to yours with 30K nodes and 800K edges) - from that exact pipeline. Took ~10 minutes read to render on my desktop:

Finally, as others have mentioned, graph partitioning might be the way to go. I really like InfoMap. It partitions the graph into a hierarchical structure that you can explore top down.

That's certainly a lot of nodes and edges, however, there is a preference called viewThreshold that should prevent the creation of a view if the number of nodes+edges is larger than the set value.  I've just tested it on a large network I have (9000 nodes and 4,000,000 edges) and it doesn't appear to be working properly.  I'll mark this as a bug.  It does, however, render, but it's not really a useful network on Cytoscape 3.7.  It will actually be usable on Cytoscape 3.8 based on early testing, however.  If you want to make your network available somewhere, I'm happy to try to see why it's not loading.
-- scooter

Actually, after further testing, the viewThreshold property works fine.  I was testing with a session file, not an import.


text:
cn.MOPS fails with 'missing value where TRUE/FALSE needed'
I am trying to use cn.MOPS to call CNVs on a set of whole exome sequencing data bam files. My script (up to the point where it fails) is:
source("https://bioconductor.org/biocLite.R")
library(cn.mops)
BAMFiles <- c('sample1.bam', 'sample2.bam', 'sample3.bam', 'sample4.bam', 'sample5.bam', 'sample6.bam', 'sample7.bam')
segments <- read.table("region.bed",sep="\t",as.is=TRUE)
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
X <- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,parallel=11)
resCNMOPS <- exomecn.mops(X)

The last command, resCNMOPS <- exomecn.mops(X) fails with:
Error in if (all(segMedianT == 0)) { :
  missing value where TRUE/FALSE needed
In addition: Warning message:
In normalizeChromosomes(X, chr = chr, normType = normType, qu = normQu,  :
  Normalization for reference sequence  chrUn_gl000228 not applicable, because of low number of segments

If I run traceback immediately afterwards, I get:
> traceback()
4: FUN(newX[, i], ...)
3: apply(sINI[chrIdx, , drop = FALSE], 2, segment, minSeg = minWidth,
       ...)
2: cn.mops(input = input, I = I, classes = classes, priorImpact = priorImpact,
       cyc = cyc, parallel = parallel, normType = normType, normQu = normQu,
       norm = norm, sizeFactor = sizeFactor, quSizeFactor = quSizeFactor,
       lowerThreshold = lowerThreshold, upperThreshold = upperThreshold,
       minWidth = minWidth, segAlgorithm = segAlgorithm, minReadCount = minReadCount,
       useMedian = useMedian, returnPosterior = returnPosterior,
       ...)
1: exomecn.mops(X)

I am using R version 3.4.4 on an Ubuntu 16.04.5 LTS machine.
My bed file only has one region which is on chrUn_gl000228:
$ grep chrUn_gl000228 region.bed
chrUn_gl000228  112604  113879

If I remove that region and rerun the same code, it works as expected. So, my question is, how can I detect this and fix it automatically? Given a collection of bam files and the bed file describing the capture regions used to sequence the samples and generate the bam files, how can I check if cn.MOPS will complain about any of the segments?
Presumably, once I have built these two objects:
segments <- read.table("region.bed",sep="\t",as.is=TRUE)
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
X <- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,parallel=11)

I should be able to then check them and see if they have the right data and, if not, remove the offending segment from the GRanges object. But how can I do this?

You start your question with biocLite, which indicates that this is not the latest version of R and Bioconductor. Which versions are you using? Have you tried with the most up to date version? Why does it complain that you have low number of segments? Could you post how many segments do you have on X?

@llrs no, I haven't used the latest, but upgrading is not an option (this is a large server). I would be happy to post that, but can you give me a command for it? My R knowledge is extremely limited. That sort of command is exactly the kind of thing I am asking about. I have no idea how to inspect and get statistics about the contents of a GRanges object.

Also, hang on, why are you saying it isn't the latest version? I got that code from the documentation of cn.MOPS. Is there a better way to load biocLite?

Well you will need to use sessionInfo() and paste here the result to see the version of R and Bioconductor. The most up to date documentation of cn.MOPS is on the bioconductor page AFAIK (if packages are not up to the standards they are kicked out of the project). Try with summary(X) and dim(X) or simply print the object and post the first lines

I said you weren't using the latest version because in the latest version the recommended way to install packages changed to BiocManager::install("cn.Mops") you need a package in CRAN to install packages of Bioconductor

@llrs yes, I added the R version (3.4.4). I had already tried summary(X) but that only prints "GRanges object with 191743 ranges and 7 metadata columns".  dim(X) gives NULL.

Let us continue this discussion in chat.

I managed to track down what was producing this error. The problem was that my bed file only had one segment for chromosome chrUn_gl000228. However, the cn.MOPS method relies on normalizing read counts across multiple segments for each chromosome. Since this chromosome only had a single segment, the exomecn.mops(X) was failing.
So, either remove the single segment from the bed file, or use a bed file that has at least 2 segments on each chromosome. I created a fake segment on chrUn_gl000228 and confirmed that the program worked as expected and didn't throw any errors.


text:
lower mapping rates in salmon v0.13 compared to previous versions
Hi there :) Thanks for the tool!
I recently updated to the new salmon (from 0.8... its been a couple years) and I noticed that my mapping percentages change dramatically between the two versions.
For example, using the default settings in v0.8, I see a mapping rate of 96.25% however upon upgrading to v0.13, these mapping rates go down to 82-84% depending on which of the parameters I play with (I've permutated  the validateMappings, incompatPrior, consensusSlack, minScoreFraction, ma, etc). Below I have a snippet of one of these tests.  I'm currently assembling the unmapped reads to see if I can determine why these reads were not mapped.
What do you think is the most likely cause for this dramatic change in performance?
Perhaps it could be the transcript index? (I used a kmer of 31 and my reads are 100 bp).
UPDATE:
I have run the indexing and quasi-mapping  using salmon v0.8, v0.10 (when validateMappings was released), v0.12, v0.13 and had mapping rates of 96, 90, 90 and 83 respectively.  So it looks like a gradual decay of mapping rates.  Or alternatively, a gradual reduction in spurious mapping?
Seems like what ever improvements from 0.12 to 0.13 changed the mapping rates significantly.  I will go ahead with both mappings (v0.12 and v0.13) to see if there are differences in the DESEQ2 outputs.
Thanks a lot!
(while I have presented values here for one set of fastq files, the same pattern is seen throughout all 72 of my samples from different experiments)
for F in $FILES ; do
  R1=${F}_L001_R1_001.fastq.gz
  R2=${F}_L001_R2_001.fastq.gz
  salmon quant \
      -p 35 \
      -i $ref_dir/$transcript_index \
      -l A \
      -1 $R1 \
      -2 $R2 \
      -o $(basename $F).VM.rFB4.IP9.quant \
      --incompatPrior 9.9999999999999995e-21 \  # default of v0.8
      --writeUnmappedNames \
      --validateMappings \
      --rangeFactorizationBins 4

done

Update salmon is now released with a --allowDovetail option in  0.13.1 and inspected via a -z dump.

Comments have been moved to chat; please do not continue the discussion here. Before posting a comment below this one, please review the purposes of comments. Comments that do not request clarification or suggest improvements usually belong as an answer, on [meta], or in [chat]. Comments continuing discussion may be removed.

[Most likely answer based on comments attached to the question]
Hi Courtney, We've looked into the data and identified the source of the different mapping rates. Specifically, the cause is discordant reads; in this case reads that are dovetailing with respect to their positions on the reference. If you run Bowtie2 with "RNA-seq" flags (i.e. disallowing discordant mappings), then you get essentially the same mapping rate as salmon here. Is allowing dovetailing reads something that should be very important in your analysis? If so, we can consider what would be necessary to add back support, but generally, they provide evidence of tenuous quality.


text:
Error in formation of correlation matrix-(corrupt matrix -- dims not not match length)
I am working on a script to create a correlation matrix and mutual_rank matrix from RNAseq gene counts for 57992 genes in 7027 samples. I have already tested this script for a smaller file (for 10,000 genes) and it is working, however, when I am running this script on my final input file, it is giving the following error:
 corrupt matrix -- dims not not match length (this error comes after gene to
 gene correlation step)

Script which i am running is:
E <- read.table("mdat_human_test.txt", header = TRUE, row.names=1)

# gene-to-gene correlation
C <- cor (t(E))
write.table(C, "human_correlation_genes_mdat.txt", quote=FALSE, sep="\t",
row.names=TRUE)

# correlation rank matrix
R <- t(apply(1-C, 1, rank))
write.table(R, "human_rank_genes_mdat.txt", quote=FALSE, sep="\t",
row.names=TRUE)

# mutual rank matrix
M <- sqrt(R * t(R))
write.table(M, "human_mutual_rank_genes_mdat.txt", quote=FALSE, sep="\t",
row.names=TRUE)

q()

What is your biological question of interest ? Why are you doing this correlations? In addition to Devon's answer, I would like to point out that if you are using all the human annotated genes you might need to filter and normalize them. You'll presumably will have 25% of "genes" which are poorly described and/or with few counts are you really interested in keeping those ?

You can't compute a correlation matrix when you have more than ~46000 rows, since a standard R matrix can have a maximum of 2^31-1 values. Have a look at packages like bigcor. Alternatively, consider if having 3+ billion correlation values is really useful.


text:
Getting data from fastq by generator
I have a task in a training that I have to read and filter the 'good' reads of big fastq files.
I downsampled, got the code working, saving in a python dictionary. But turns out the original files are huge and I rewrite the code to give a generator. It did work for the down-sampled sample. But I was wondering if it's a good idea to get out all the data and filtering in a dictionary. Does anybody here a better idea? I am asking because I am learning python on my own and I am not sure if it makes sense.
I got some ideas from a code in Biostar:
import sys
import gzip

filename = sys.argv[1]

def parsing_fastq_files(filename):

    with gzip.open(filename, "rb") as infile:
        count_lines = 0
            for line in infile:
            line = line.decode()
            if count_lines % 4 == 0:
                ids = line[1:].strip()
                yield ids
            if count_lines == 1:
                reads = line.rstrip()
                yield reads
        count_lines += 1

total_reads = parsing_fastq_files(filename)
print(next(total_reads))
print(next(total_reads))

I now need to figure out to get the data filtered by using 'if value.endswith('expression'):' but if I use a dict for example, but that's my doubt because of the number of keys and values.

You want lines where the header endswith expression correct?

Have a function that returns a single-end read (or pair) and another that does any relevant filtering on that. Wrap both of them in another function that checks if you've hit the end of the file. That's both simple, will handle any size of file, and allows single and paired-end reads.

This is more of a comment than an answer right Devon?

@d_kennetz Not really, it answers what the proper scalable method is. One could include the actual code, but I have faith in the poster's ability given a prose description.

I do not have an example of fastqs with 'expression' as the last string in the header, so instead I will filter based off of a specific barcode used for demultiplexing. I have a gzipped fastq of 250 reads in size. Based on Devon's answer, I have first generated a fastq parser which stores only (header, seq) as a list of tuples, and then looks inside the list of tuples for a specific index. Some of the indexes vary by 1 nucleotide so those will be left out. I can prove this by printing the length of the 2 lists.
Here is the code:
import sys
import gzip

fastq = sys.argv[1]

def parser(filename):
    """
    Generate a list of tuples (header, read)
    """
    fastq_parsed = []
    try:

        with gzip.open(filename) as fq:
            header = next(fq)
            read = next(fq)
            fastq_parsed.append((header[1:-1], read[:-1])) # we don't want new-line chars or @ sign in reads
            while True:
                next(fq) # skip + line
                next(fq) # skip qscore
                header = next(fq) # store next read header
                read = next(fq) # store next read seq
                fastq_parsed.append((header[1:-1], read[:-1]))
    except:
        StopIteration # when we run out of reads, stop
    return fastq_parsed

def expression_finder(read_list):
    """
    Find headers that end with expression assuming the list format generated above.
    """
    filtered_reads = []
    for header, read in read_list:
        if (header.endswith(b'ATCACGAT')): #gzip files in bytes format, so we pass bytes object.
            filtered_reads.append((header.decode("utf-8"), read.decode("utf-8"))) # list items decoded to utf-8
    return filtered_reads

read_list = parser(fastq)
print(len(read_list))
subset_reads = expression_finder(read_list)
print(len(subset_reads))

Which returns:
(base) [dkennetz@nodecn203  test]$ python3.6 fastq_parser.py subset.fastq.gz
250
173

We can see that of the 250 reads, 173 have the correct index. I have hardcoded the filter in the second function, but you could also pass it as an argument to the function, like:
def expression_finder(read_list, filter):
    ...
        if header.endswith(filter):
            ...

subset_reads = expression_finder(read_list, b'expression')

which would also work. But include the b'' in your string syntax because the first function still passes a bytes object.

Hey d looks cool but if you .decode() the header and the read when you append? Thanks for share your ideas. paulo

I did decode when I append but you certainly don’t have to do so! You can leave them as bytes if you’d like by simply passing .append((header, seq)). This will leave each as bytes objects.


text:
Assembly by stringtie
I run this cmd
./stringtie G1_sorted.bam -B -o G1.gtf -G Triticum_aestivum.IWGSC.42.gtf -p 4  -C G1.refs.gtf -A G1.abund.tab

Error is:
-WARNING: no reference transcripts were found for the genomic sequences where reads were mapped!
Please make sure the -G annotation file uses the same naming convention for the genome sequences.

How can I resolve it?
Update:
I have downloaded both (GTF and Refernce Genome)files from Ensemble. Can you confirm me that whether theu have same chromosomes naming? So that i can perform my mapping  task of GTF File with Reference Genome which usually i remain fail to do due to difference in chromosomes naming.                                                           Here The Reference genome file is shown as:
 >1A dna:chromosome chromosome:IWGSC:1A:1:594102056:1 REF
CTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCCTAACCCTAAACCCTAACCAAAAC

Where as The annotation file (gtf file ) is as
    #!genome-build IWGSC
#!genome-version IWGSC
#!genome-date 2018-07
#!genome-build-accession GCA_900519105.1
3B      IWGSC   gene    212892  214491  .       -       .       gene_id "TraesCS3B02G000100"; gene_source "IWGSC"; gene_biotype "protein_coding";
3B      IWGSC   transcript      212892  214491  .       -       .       gene_id "TraesCS3B02G000100"; transcript_id "TraesCS3B02G000100.1"; gene_source "IWGSC"; gene_biotype "protein_coding"; transcript_source "IWGSC"; transcript_biotype "protein_coding";

Can it be confirmend that these both files are same ? If not same with chromosome naming convention , then please suggest me how can I change it ?

This looks like your chromosome names in the BAM don't match the GTFs you're using. Check that first.

the reference genome download from the NCBI.and gtf  file from ensembl. your mean is that i have to use the gtf file from the NCBI, then this error will be removed?

Welcome to the site @fatima. Yes ENSEMBL uses genes identifiers starting with ENSGXXXXX.XX while NCBI uses gene names with just numbers and they don't map to the same positions. They aren't exacts equivalent

hi IIrs, i have use the reference genome and gtf file from ensembl.encounter same problem.how i can mofify the file with same name?

hi @fatima (to send me a notification use @ followed by my username) Could you please explain how do you want to modify these files to use the same name? rereading again your question I find it a bit unclear what do you want to achieve and what do you expect us to do for you.

This seems to be a common problem. As @llrs pointed out, Ensembl, UCSC(genome browser and associated data) and NCBI records are not directly interchangeable. The ambiguity refers to annotation as well as the reference genome. The reference genome is identical, except for the chromosome naming convention, meaning a reference from Ensembl might not contain CHR prefix, whereas it will be present in one downloaded from UCSC.
To prevent such conflicts in downstream tools, it is good to stick to the reference and annotation from the same repository throughout the analysis and map the genes to gene names or whatever towards the end of the investigation.
You can download all the associated data for T. aestivum which I beleive is the organism you want to study from Ensembl plants FTP download.

Hi Siddharth, Thank you, Now I have downloaded both files from Ensembl. Can you please tell me these are with the same naming?                                                                                I have open the reference file.The Reference genome file is as                                                >1A dna:chromosome chromosome:IWGSC:1A:1:594102056:1 REF
CTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCCTAACCCTAAACCCTAACCAAAAC
The annotation file (gtf file ) is as


text:
bedtools unable to open file despite being tab-delimited
I got the following error when attempting to intersect two files. The file below is supposed to be the -a file.
Error: Unable to open file compare_me_to.txt. Exiting.

When I attempted to make it tab-delimited through Ubuntu (see code below)
perl -p -i -e 's/ /\t/g' compare_me_to.txt

I got the following error:
Can't open compare_me_to: Permission denied.

Additionally, the ownership of the folder the files are in is of my username.
I can't tell if this is a bedtools problem or a Ubuntu problem or something else.

Just a quick check since you said you have ownership of the file... The file isn't open anywhere else (like in another shell or on desktop) when you try to convert right?

I suspect the reason is the "-i". The following will work,
perl -pe 's/  */\t/g' compare_me_to.txt > compare_me_to.txt2

-i will want to write onto the file and create a backup. I suspect you'll have spotted chmod +w compare_me_to.txt.
Just make sure that a double space is inserted before the * character, otherwise it will go wrong. Copying and pasting will avoid this.

This is a clearer solution because it avoids missing a space character,
perl -e 's/[ ][ ]*/\t/g' compare_me_to.txt > compare_me_to.txt2

and the solution you used .. (I did on OS X), what I'd recommend is,
perl -p -e 's/[ ][ ]*/\t/g' compare_me_to.txt # to see if it works by printing on screen without altering the file

perl -i -e 's/[ ][ ]*/\t/g' compare_me_to.txt

The  -i flag will alter the original file and store a backup without requiring renaming the file.
What's it doing?
The regex [ ][ ]* structure is saying look a space and if one or more spaces occur adjacent to each other replace everything with a single tab. [ ]* means a "space" that is present or absent and if present it can be present an infinite number of times. If the [ ]* is by itself the command risks saying replace '' (zip) or 'space' with a tab, so tabs can appear at every '' (zilch), so tabs can appear everywhere. The "*" symbol in regex work differently to the "*" wildcard character in shell.

The following error appears: Can't open compare_me_to.txt: Permission denied.

I was able to give myself permission (https://askubuntu.com/questions/466605/cannot-open-output-file-permission-denied), but the file still can't be opened when intersecting. I ran the line of code you wrote; when I opened the resulting file, every character was separated by a tab.

Oh I know, you put a single space into the term ... its a double space followed by a *. so its ... 's/  */\t/g' ... if you only put a single space into the first / */ then it will result in tabs going everywhere. I can explain later.

chmod 700 file_name is another way in, I'm surprised chmod +w file_name didn't work if you were the file owner. I assumed you knew about chown and chmod


text:
What is the purpose of 'Folder locations' in MaxQuant?
MaxQuant is a tool for proteomics quantification. In the settings under Global parameters / Folder locations are three folder locations that are possible to be changed. Thought the manual does not even mention these options. Is there any benefit in changing these locations? Or can it be bad to do that? For example, if I run multiple MaxQuant instances at the same time doing different searches?
What is the meaning of the Fixed andromeda index folder, Temporary folder and Combined folder location in MaxQuant? What is it used for and when would you (not) change it?

These folders are usually generated in a default directory (where the raw files are). It can however be beneficial to set the paths to these folders manually.

Temporary folder

that is obviously where temporary files go, so having an SSD or other fast drive will benefit the I/O operations. Set that path to a fast storage partition.

Combined folder

all the output, including txt directory will go here. Setting it to fast storage might speed up writing the txt files. If you set it to the same path while running multiple instances of MQ it will obviously override output.

Andromeda index

this is where index files from Andromeda search are. By default these are inside a combined/search folder. Since these index files can get big, once generated by MQ they can be reused. Set this path if you want to reuse previously generated index files.


text:
maf file format: Variant_Classification value for non frameshift substitution
I have a vcf containing SNV and MNV (short indells or substitutions), I have annotated that file with annovar and in the resulting txt file I have this row:
chr7    140453136   140453137     AC   TT  exonic   BRAF    .     nonframeshift substitution              BRAF:NM_004333:exon15:c.1798_1799AA:p.V600K

If I convert the txt file in maf using annovarToMaf function of maftools bioconductor package I obtain NA as value for Variant_Classification field. In other words, the field is included in the original vcf file but it is converted into NA when saved as maf with maftools.
What is the correct value of the Variant_Classification field in a .maf file for non-frameshift substitutions like BRAF V600K?
I tried vcf2maf and that tool put Missense_Mutation as a value for Variant_Classification of BRAF V600K mutations.
I tried to download the .maf file of TCGA melanoma (that should contain BRAF V600K mutations) but the access to that file is restricted.

Ok, try to make a reproducible example as it might be a bug on the maftools package. However I would first post this reproducible example here to see if there is a solution without involving changes on the package


text:
HISAT2 - are all multi mapped reads reported by default?
While going through the HISAT2 manual, one can see this:

Two alignments for the same individual read are "distinct" if they map
  the same read to different places. Specifically, we say that two
  alignments are distinct if there are no alignment positions where a
  particular read offset is aligned opposite a particular reference
  offset in both alignments with the same orientation. E.g. if the first
  alignment is in the forward orientation and aligns the read character
  at read offset 10 to the reference character at chromosome 3, offset
  3,445,245, and the second alignment is also in the forward orientation
  and also aligns the read character at read offset 10 to the reference
  character at chromosome 3, offset 3,445,245, they are not distinct
  alignments.
Two alignments for the same pair are distinct if either the mate 1s in
  the two paired-end alignments are distinct or the mate 2s in the two
  alignments are distinct or both.
Default mode: search for one or more alignments, report each

and then there is -k parameter (analogous to bowtie2 -K parameter)
which defines:

HISAT2 searches for up to N distinct, primary alignments for each
  read, where N equals the integer specified with the -k parameter.
  Primary alignments mean alignments whose alignment score is equal or
  higher than any other alignments. It is possible that multiple
  distinct alignments have the same score. That is, if -k 2 is
  specified, HISAT2 will search for at most 2 distinct alignments. The
  alignment score for a paired-end alignment equals the sum of the
  alignment scores of the individual mates. Each reported read or pair
  alignment beyond the first has the SAM 'secondary' bit (which equals
  256) set in its FLAGS field. See the SAM specification for details.

This is a bit confusing, as the first quote means that HISAT2 finds all distinct valid alignment and reports them by default, and the second quote says that when -k option is used, it will only report a subset of alignments which it would have reported by default. Am I correct in this inference or does one needs to use -k option to enable "multi-mapping"?

You said:

...the first quote means that HISAT2 finds all distinct valid alignment and reports them by default.

Not exactly. The first section you quoted doesn't say it will find and report all distinct alignments, but rather that it would find and report one or more distinct alignments. It's a subtle difference, but an important one.
That section describes how distinct mappings are defined, not how many are reported. The second section you quoted, as well as the documentation for the -k parameter, are all consistent with this interpretation.

-k <int>        It searches for at most <int> distinct, primary alignments for
                each read. Primary alignments mean alignments whose alignment
                score is equal or higher than any other alignments. The search
                terminates when it can't find more distinct valid alignments, or
                when it finds <int>, whichever happens first. The alignment
                score for a paired-end alignment equals the sum of the alignment
                scores of the individual mates. Each reported read or pair
                alignment beyond the first has the SAM 'secondary' bit (which
                equals 256) set in its FLAGS field. For reads that have more
                than <int> distinct, valid alignments, hisat2 does not guarantee
                that the <int> alignments reported are the best possible in
                terms of alignment score. Default: 5 (HFM) or 10 (HGFM)

                Note: HISAT2 is not designed with large values for -k in mind,
                and when aligning reads to long, repetitive genomes large -k can
                be very, very slow.

By default, 5 or 10 distinct mappings are reported (depending on the mode), but you can control this with the -k parameter.

Aha, thanks for clearing it up! Although I'm still a bit confused to why HISAT2 still has --secondary if setting up -k instructs it to search for multiple alignments n alignments regardless.

The -k parameter corresponds only to primary alignments: cases where distinct alignments have an identical high score. If an alignment has a less-than-optimal score, it is a secondary alignment.


text:
Is there a resource where I can donwload E. coli assemblies grouped by their pathotypes?
I will preface this by saying I'm not a microbiologist, so I apologise for any haziness and I'll be happy to expand.
I would like to download the assemblies of different E. coli genomes and group them by their pathotypes (e.g. EPEC, ETEC, EIEC).
I have used enterobase previously to download salmonella genomes broadly by their subspecies, and it does have a set of E. coli genomes. However, they are mainly annotated by their serotype, and not with any broad nomination of pathogenesis. I think I would be able to use the API to grab these serotypes and map them to a type of pathogenicity but that would probably be quite messy.
Does such a database exist? Most I can find are broadly annotated or a are very detailed, but don't allow access to the raw annotation.
Thanks in advance.

If ENA archive contains such information, then you can access it programatically through the EBI api.
Also you can access the NCBI database programatically through their eutils.
In both cases you can access annotations, though not sure if they contain what you are looking for in all the cases.

Thanks for that! I previously tried the NCBI and it turned out to be quite messy. The ENA is a step in the right direction as I can extract the strain information quite nicely. Not quite what I was after but definitely helpful, thanks.


text:
Is there a risk in inference of structural variations using mate-pairs?
Inference of structural variants from short read data is always sneaky (for a question about pair-end reads, see this). Recently I bumped into a paper where they used mate-pairs for inference of structural variants and it made me wonder if the approach might be problematic or not.
Mate-pairs have a known problem of contamination by pair-end reads and there is even a tool for calling structural variations considering this problem called SVachara. However, authors have used methods Lumpy and Delly2. Have they chosen an appropriate approach?


text:
Picard validation error regarding bin field of BAM file
I am sure I set the path right but whatever I am trying the command not working
Any help please?
[fi1d18@cyan02 fi1d18]$ picard MarkDuplicates \
    I=DNA.dupmarked.bam \
    O=DNA.dupmarked1.bam \
    M= marked-dup-metrics.txt
[Thu Mar 07 17:33:42 GMT 2019] picard.sam.markduplicates.MarkDuplicates INPUT=[/temp/hgig/fi1d18/1631_WTSI-OESO_005_a_DNA/mapped_sample/HUMAN_1000Genomes_hs37d5_genomic_WTSI-OESO_005_a_DNA.dupmarked.bam] OUTPUT=/temp/hgig/fi1d18/1631_WTSI-OESO_005_a_DNA/mapped_sample/HUMAN_1000Genomes_hs37d5_genomic_WTSI-OESO_005_a_DNA.dupmarked1.bam METRICS_FILE=marked-dup-metrics.txt MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP=50000 MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 SORTING_COLLECTION_SIZE_RATIO=0.25 REMOVE_SEQUENCING_DUPLICATES=false TAGGING_POLICY=DontTag REMOVE_DUPLICATES=false ASSUME_SORTED=false DUPLICATE_SCORING_STRATEGY=SUM_OF_BASE_QUALITIES PROGRAM_RECORD_ID=MarkDuplicates PROGRAM_GROUP_NAME=MarkDuplicates READ_NAME_REGEX=<optimized capture of last three ':' separated fields as numeric values> OPTICAL_DUPLICATE_PIXEL_DISTANCE=100 VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_INDEX=false CREATE_MD5_FILE=false GA4GH_CLIENT_SECRETS=client_secrets.json
[Thu Mar 07 17:33:42 GMT 2019] Executing as fi1d18@cyan02 on Linux 2.6.32-754.el6.x86_64 amd64; Java HotSpot(TM) 64-Bit Server VM 1.8.0_51-b16; Picard version: 2.8.3-SNAPSHOT
INFO 2019-03-07 17:33:42 MarkDuplicates Start of doWork freeMemory: 2012347496; totalMemory: 2027945984; maxMemory: 3817865216
INFO 2019-03-07 17:33:42 MarkDuplicates Reading input file and constructing read end information.
INFO 2019-03-07 17:33:42 MarkDuplicates Will retain up to 14684096 data points before spilling to disk.
WARNING: BAM index file /temp/hgig/fi1d18/1631_WTSI-OESO_005_a_DNA/mapped_sample/HUMAN_1000Genomes_hs37d5_genomic_WTSI-OESO_005_a_DNA.dupmarked.bam.bai is older than BAM /temp/hgig/fi1d18/1631_WTSI-OESO_005_a_DNA/mapped_sample/HUMAN_1000Genomes_hs37d5_genomic_WTSI-OESO_005_a_DNA.dupmarked.bam
[Thu Mar 07 17:33:42 GMT 2019] picard.sam.markduplicates.MarkDuplicates done. Elapsed time: 0.01 minutes.
Runtime.totalMemory()=2027945984
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
Exception in thread "main" htsjdk.samtools.SAMFormatException: SAM validation error: ERROR: Record 3752, Read name HX3_22030:3:2114:23155:23319, bin field of BAM record does not equal value computed based on alignment start and end, and length of sequence to which read is aligned
at htsjdk.samtools.SAMUtils.processValidationErrors(SAMUtils.java:448)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.advance(BAMFileReader.java:665)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.next(BAMFileReader.java:650)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.next(BAMFileReader.java:620)
at htsjdk.samtools.SamReader$AssertingIterator.next(SamReader.java:569)
at htsjdk.samtools.SamReader$AssertingIterator.next(SamReader.java:543)
at picard.sam.markduplicates.MarkDuplicates.buildSortedReadEndLists(MarkDuplicates.java:438)
at picard.sam.markduplicates.MarkDuplicates.doWork(MarkDuplicates.java:222)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:205)
at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:94)
at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:104)
[fi1d18@cyan02 fi1d18]$

The error message is:

Exception in thread "main" htsjdk.samtools.SAMFormatException: validation error: ERROR: Record 3752, Read name HX3_22030:3:2114:23155:23319, bin field of BAM record does not equal value computed based on alignment start and end, and length of sequence to which read is aligned
EDITED

[fi1d18@cyan02 sambamba]$ samtools view -H /temp/hgig/fi1d18/1631_WTSI-OESO_005_a_DNA/mapped_sample/HUMAN_1000Genomes_hs37d5_genomic_WTSI-OESO_005_a_DNA.dupmarked.bam
[W::bam_hdr_read] EOF marker is absent. The input is probably truncated.
@HD     VN:1.5  SO:coordinate
@SQ     SN:1    LN:249250621    AS:NCBI37       M5:1b22b98cdeb4a9304cb5d48026a85128     SP:Human  UR:/lustre/scratch117/core/sciops_repository/references/Homo_sapiens/1000Genomes_hs37d5/all/fasta/hs37d5.fa
@SQ     SN:2    LN:243199373    AS:NCBI37       M5:a0d9851da00400dec1098a9255ac712e     SP:Human  UR:/lustre/scratch117/core/sciops_repository/references/Homo_sapiens/1000Genomes_h

The java executable is not on your PATH.

Thank you but when I loaded likely related module I got another error in edited post

Can you post what that line of the bam looks like?  Did you try remaking the bam index, since the software complained about it?

Have you considered just using sambamba to mark duplicates? It tends to be faster and less picky. There were a couple old versions of samtools that were messing up the bin assigned to some records, so make sure you use the most recent version in the future if you're not already doing so.

Sorry I pasted a few lines of my bam in my main post please consider

Sorry @Devon sambamba says    [fi1d18@cyan02 sambamba]$ make sambamba-ldmd2-64
make: *** No rule to make target `sambamba-ldmd2-64'.  Stop.

Why would you try to compile programs rather than using conda or downloading the binary?

Please read error messages thoroughly, the cause is written right in what you posted.

[fi1d18@cyan03 ~]$ samtools quickcheck file.bam
[fi1d18@cyan03 ~]$

Let us continue this discussion in chat.

We've seen this error before on the GATK forum; if you google "bin field of BAM record does not equal value computed based on alignment start and end" you will find this forum thread:
https://gatkforums.broadinstitute.org/gatk/discussion/4290/sam-bin-field-error-for-the-gatk-run
This part of the last post should help you solve this:

The only way I can think of to fix this is to fix the BAM. In order to
  do that, you would need to convert to SAM format, which does not have
  a field for indexing bin, and then back to BAM. That would be the
  safest thing, followed by recreating the BAI. So you can do this with
  GATK tools using the following commands.
gatk SamFormatConverter - to convert Bam to Sam

gatk SamFormatConverter - to convert Sam back to Bam

gatk BuildBamIndex - to create the Bam index file.

Also, please make sure your Picardtools and GATK versions are up to
  date to the most recent versions.


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:
Geneland R + Fortran error?
I am helping an undergrad with their project and their thesis advisor has requested they use GENELAND for their analysis. I haven't worked on this sort of analysis, so all I have done is install it in their R, hand over the manual, and help troubleshoot any errors they can't resolve.
We've gotten the analysis through the MCMC stage and are getting an error with the PostProcessingChain command.
Their MCMC command is:
MCMC(geno.hap=geno1, spatial=FALSE, nb.nuclei.max=3,
    varnpop=TRUE, freq.model="Uncorrelated",
    filter.null.alleles=FALSE, nit=100000,
    thinning=100, npopmin=1, npopinit=3,
    npopmax=10, rate.max=3, delta.coord=0,
    path.mcmc="/Users/myname/localdoc/GENELAND")

Their PostProcessChain command is:
PostProcessChain(coordinates=gene1coord.utm,
    path.mcmc ="/Users/myname/localdoc/GENELAND/",
    nxdom=100, nydom=100, burnin=200)

This is exactly the same command as in the manual, and gene1coord.utm is an R object with the coordinates transformed into UTM coordinates. The error we receive is:

([1] "Calling Fortran function postprocesschain2"

At line 12334 of file main.f (unit = 10)

Fortran runtime error: File already opened in another unit)

I have tried fix(PostProcessChain) to see if there's something I can adjust and this is what I see for the section that shows up:

So far, here is the list of things I have tried:

Reinstall X11
Reinstall gfortran
make sure that R is 3.4.x as the GENELAND package suggests
Tried different parameters on MCMC just in case it was user error
Tried different parameters on PostProcessChain, but they're the same as the manual really
Tried different kinds of coordinates
Tried a spatial analysis instead of a non-spatial
Googled the error
Restarted the computer

Alternatively, is there another program that does similar analysis to geneland? I'm not sure if the advisor recommended it because it's the only thing that does what it does or because that's just what they used back in the day.
Edit: There were 3 blank files in the folder titled %FO, P, and x. When I deleted them the "file open" error went away but it still hangs up on [ 1 ] Calling Fortran function postprocesschain2 and has a fatal error that crashes R.

I don't suppose the file is currently being accessed by anything else, right? Just to be 100% sure, have you tried rebooting (in case something weird happened and you have a stale open file handle somewhere)? Also, please avoid posting images of text. That makes it considerably harder for us to copy the error and look for it.

@terdon Thanks, no one/nothing else could be accessing the files but R. There isn't even a file called "main.f" or any files with 12000 lines in the folder the MCMC output goes to, and PostProcessChain accesses.
I tried a restart and it changed nothing.


text:
Does UCSC definition for telomere making sense?
UCSC provides a database of telomere regions for each chromosome:

ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/gap.txt.gz

It reported the telomere section in chr21 as:
grep chr21 gap.txt | grep telomere
585 chr21   0   10000   1   N   10000   telomere    no
952 chr21   48119895    48129895    515 N   10000   telomere    no

Telomere is the first 10000 bases in chr21?? But ... the region are all "N" in the reference genome!
How does UCSC confidently give out telomere regions if they don't even know anything about the sequences?

Have you read any of the literature on this?

It seems so, yes.
UCSC uses a 10Kb definition and representation for telomeres in human. In the table you shared this is consistent across all chromosomes, not only chr21. This 10kb representation is an approximation to the mean value observed in newborns and is supported by evidence like the one linked below.
Taken from Wikipedia:
A telomere is a region of repetitive nucleotide sequences at each end of a chromosome, which protects the end of the chromosome from deterioration or from fusion with neighboring chromosomes. For vertebrates, the sequence of nucleotides in telomeres is AGGGTT, with the complementary DNA strand being TCCCAA, with a single-stranded TTAGGG overhang. This sequence of TTAGGG is repeated approximately 2,500 times in humans. In humans, average telomere length declines from about 11 kilobases at birth[ to less than 4 kilobases in old age, with the average rate of decline being greater in men than in women.
2,500 x 6 (TCCCAA) = 15,000
(approx. 10K)
Two interesting additional facts:
The minimum telomere length needed to ensure human telomere protective stability in white blood cells is 3.81 kb
Mouse has a 100Kb telomere representation in USCS


text:
How to interpret fish plot in R
I am working on subclonal evolution of multiple myeloma. After literature survey, I found people use Fish plot to represent tumor evolution. The example of fish plot can be found here. But I am not able to interpret this plot. How this plot represent tumor evolution and subclones in tumor. Can anybody explain with example. Thanks.

Have you read the paper cited in your link? Doesn't that explain them? We can't help you understand a plot if you don't even tell us what is being plotted, what the y-axis represents.

@terdon : This is what I did as first step. But in their paper, they have not given any such explanation or interpretation of Fish plot (implementation info given only) . Paper can be found at (https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-016-3195-z). I was hoping to get some useful information in paper but not.

The population (cells) is going through a bottleneck, whether that is a population size, aneuploidy, dysregulation I've no idea, the y-axis is not labelled. If it's not an oncovirus I'm stumped. Oh it's longitudinal study BTW (obviously)

Author of the package here, stumbled across this today.  The existing explanation is close, but not entirely correct.
Our paper describes it like this: "We first created a fish plot of an AML patient with a chemotherapy-induced bottleneck that eliminated one subclone, while another survived and drove the relapse".
In other words, the grey represents a tumor, starting from a single cell (point at far left). The tumor grows and expands, acquiring additional mutations (red) that increase the relative growth rate of a subclone.  When sampled at day 0 (when the patient showed up with cancer), the red subclone made up about a third of the tumor cells.  Chemotherapy was applied, shrinking the tumor dramatically and completely eliminating the red subclone (patient went into remission).
Unfortunately, remission doesn't always mean "cured", as this tumor persisted, and picked up more mutations (orange subclone), and then within that subclone, picked up even more (yellow subclone).  Each gave those cells sufficient growth advantage that at relapse, essentially all of the tumor cells were comprised of the orange subclone, with a subset of cells that contained both the orange and yellow mutations.
There are also subtleties about whether the orange and yellow mutations pre-existed at very low levels (below the sensitivity of the sequencing that was performed), or were actually brand new mutations, but that's outside the scope of this explanation!

Each color is a clone (set of mutations) or an individual mutation (depends on how the input data was generated). This is a cancer with the bottleneck likely being induced by a treatment (e.g. UV, chemo, drugs, etc...)
Based on the vignette this is likely an AML (leukemia).
It looks like grey are healthy cells. At day 0 a specific clone (red) existed at a low frequency (~20-40%). The patient (or cells) was then exposed to a treatment and on day 150 another measurement was taken. In this case two novel clones (subclones of red) have taken over the population. They seem to co-exist or are co-dominate clones with orange being 50-60% of the population and yellow being 30-40%.
The example says "Simple plot with multiple subclones" so: the orange and yellow populations have additional (likely mutually exclusive) mutations that were not present in the red clone - but biologically also have all the mutations in the red clone. So my first assumption was right - each color is a proxy for a set of mutations.
You could have generated the mutation data (i.e. allele frequencies) from a wide variety of sources (scRNA-Seq, exome-sequencing, PCR genotyping, etc...)

Thanks for this simplistic explanation. I want to know  what does that narrow gray color line, width of each sub-region represent in the above figure (between day 0 and day 150). In the same example (https://www.rdocumentation.org/packages/fishplot/versions/0.4)  the other graph with More complex clonal structure, in day 403, blue region continue till day 403 but width reduces and green/red region starts very  late in day 403 region. Does that mean that sub-clone represented by blue region still exists in tumor on day 403 and green/red sub-clones  develop on last time span?

Similarly, in figure Multiple independent clones, many timepoints, your explaination is not perfectly applicable. One sub-clone is coming separately from healthy cells. But as per your explanation, all the sub-clones should lie in between healthy cell region (even for multiple clones).

Its easy, the graph infers treatment was counter-productive. There is an initial clonal expansion of a cancerous cell (red), within healthy tissue (grey). The cancerous cells appears to be controlled by treatment, but healthy tissue is also destroyed (bottleneck). Two new clonal expansions emerge in place of the red/grey cells post-treatment and rapidly become dominant. Thus treatment has lead to a more aggressive cancer, i.e. counter-productive. In population genetics its a "selective sweep", orange/yellows clones are likely present in the originally but couldn't become dominant.


text:
Including Picard tool in galaxy
I am running the Picard tool CollectSequencingArtifactMetrics in galaxy (locally in my pc).
The output files are saved into disk and not passed to stdout and (I guess) they can't be passed as output from picad to galaxy .
In the xml file I am including:
<outputs>
    <data name='preSummary' format="tabular"  from_work_dir="pref.pre_adapter_summary_metrics"/>
</outputs>

When I click on View details, in "stderr", I see:

[Mon Mar 11 13:43:07 EDT 2019]
  picard.analysis.artifacts.CollectSequencingArtifactMetrics done.
  Elapsed time: 0.19 minutes. Runtime.totalMemory()=3323461632

1 - It looks like it worked but shouldn't I see the highlighted information in "stdout" rather than in the "stderr"?
2- If it worked, How can I get my output files? I tracked the working directory printing pwd from the <command> but I guess all files are deleted immediately once the picard process has finished.

Any idea of what am I doing wrong?
What basic knowledge am I missing here?
Thank you!

Picard, like most tools, prints logs to stderr.
It should be in the Galaxy history. Are you using a version of picard from the toolshed (if so, which one)? You might set cleanup_job to never in the config, which will prevent galaxy from cleaning the working directory. This will enable seeing what the job actually produced.

Most likely there's more in the log file that you're not showing is us and picard itself never ran to completion. In lieu of that information, what I wrote above is usually sufficient to debug things.

Picard was running (exit status was 0 when I set to 1 if there was an error) The log was the same as when I run the command in the terminal. I just excluded the prefix "./" in naming the output and it worked (I guess galaxy add a path to the filename with is not "./" ?).
Thanks!


text:
How can I specify an email address in the SLURM snakemake json configuration file?
I am working on a SLURM-based cluster. I am currently using the following slurm json file:
{
    "__default__" :
    {
        "A" : "overall",
        "time" : "96:00:00",
        "nodes": 1,
        "ntasks": 1,
        "cpus" : 24,
        "p" : "LONG",
        "mem": 16G,
        "output": "snakemake%A.out",
        "error": "snakemake%A.err"
    }
}

along with the following command:
snakemake --jobs 200 --use-conda --printshellcmds  --cluster-config slurm.json --cluster "sbatch -A {cluster.A} -p {cluster.p} -t {cluster.time} --output {cluster.output} --error {cluster.error} --ntasks {cluster.ntasks} --cpus-per-task {cluster.cpus} --mem {cluster.mem}"

I would like to be able to specify my email and receive an email when a job finishes. I have tried adding a line in the json file with "email": my@email.com and an option to the sbatch command --email {cluster.email} but that did not work.
Anyone knows how to do it?

After looking at some headers of files (#SBATCH --mail-user=email@comI have used --mail-user instead. I'll let you know if I receive an email.

I'd be surprised if --mail-user didn't end up working. Note also --mail-type.

so --mail-user on its own didn't work but by adding the --mail-type, it worked!

As noted by the comments and the SLURM manual (not so much bioinformatics related?)
You need the following terms in your batch/queuing command:
--mail-type=ALL
--mail-user=you@email.com

I used "ALL" which will send you all reports. If you just want the finish message then you would use "END" and "FAIL" like such:
--mail-type=END,FAIL

Here are the specifics from the manual if you wish to make a different selection:
--mail-type=

Notify user by email when certain event types occur. Valid type values
  are NONE, BEGIN, END, FAIL, REQUEUE, ALL (equivalent to BEGIN, END,
  FAIL, REQUEUE, and STAGE_OUT), STAGE_OUT (burst buffer stage out and
  teardown completed), TIME_LIMIT, TIME_LIMIT_90 (reached 90 percent of
  time limit), TIME_LIMIT_80 (reached 80 percent of time limit),
  TIME_LIMIT_50 (reached 50 percent of time limit) and ARRAY_TASKS (send
  emails for each array task). Multiple type values may be specified in
  a comma separated list. The user to be notified is indicated with
  --mail-user. Unless the ARRAY_TASKS option is specified, mail notifications on job BEGIN, END and FAIL apply to a job array as a
  whole rather than generating individual email messages for each task
  in the job array.

--mail-user=

User to receive email notification of state changes as defined by
  --mail-type. The default value is the submitting user.


text:
Tools for quality trimming at 5 prime?
I'm looking for a tool that is capable to do quality trimming at 5' end and it is configurable. For example, I can choose the quality trheshold, the read length, etc. Any recommendations?
I'm currently using sickle for quality trimming. However, one of my samples have very low quality at the beginning as well and sickle doesn't seem to be able to deal well with it.

Hi there, looking at README of sickle "You can also change the length and quality thresholds for trimming, as well as disable 5'-trimming and enable truncation of sequences with Ns." which sort of suggests that it actually trims 5'. I feel like I am missing something.

Yes, I've tried that, but I've got the feeling it's too stringent for this particular sample.

Ah, now I get it. So you search for a tool that is configurable for quality trimming of 5' end, right? I would recommend [edit] your question and explain this.

You can have a look at cutadapt. It is capable of quality trimming for both ends as you can read here.

I've just tried this, however I still see the 5' quality still very low. I run it like this:
cutadapt -q 30,30 -o output.fastq input.fastq
Is there any thing else I'm missing?

Could you please give some quality summary information of your reads? Before and after the trimming?

I never spent too much time on choosing my trimming software, therefore I might have missed some jewels. I use trimmomatic when I need versatility and I am entirely sure that it trims reads on both ends. A more convenient option is skewer, they do perform quality filtering, but I could not find an explicit confirmation that it works on 5' too.

Out of all of the major trimming tools available and widely used (trimmomatic, cutadapt/trimGalore (trimGalore is built on top of cutadapt), fastp), I actually instead prefer bbduk which is part of the bbmap suite:
https://github.com/BioInfoTools/BBMap
https://github.com/BioInfoTools/BBMap/blob/master/sh/bbduk.sh
It's much more "fine tuneable" than the other tools and runs in much less time using much less memory (as long as you set the memory limit as it's written mostly in java!). You can trim both ends or use a sliding window approach similar to what you see in Trimmomatic. Even better you can limit the number of bp in that it will trim from either left/right (5'/3'). It's by far the best all-around trimming tool for any type of sequencing dataset ({RNA, ChIP, PRO, GRO, ATAC, DNA}-seq etc...) in my opinion.


text:
Running htseq-count over BAM files
I am trying to derive an expression matrix from BAM files using htseq-count. These are bulk RNASeq BAM's by the way.
I have read the htseq-count documentation as well as samtools and figured that the following command should work:
samtools sort -on file.bam | samtools view -h - | htseq-count --type=exon --idattr=gene_id --additional-attr=gene_name - Homo_sapiens.GRCh38.77.gtf > counts.txt

The .gtf file here is from ensembl.
I have tried many variations of the above command but have failed to generate a count matrix. The error message states:
"0 SAM alignments processed"
Here are the first few lines of my .gtf file:
#!genome-build GRCh37.p13
#!genome-version GRCh37
#!genome-date 2009-02
#!genome-build-accession NCBI:GCA_000001405.14
#!genebuild-last-updated 2013-09
1       pseudogene      gene    11869   14412   .       +       .       gene_id "ENSG00000223972"; gene_name "DDX11L1";
gene_source "ensembl_havana"; gene_biotype "pseudogene";
1       processed_transcript    transcript      11869   14409   .       +       .       gene_id "ENSG00000223972"; trans
cript_id "ENST00000456328"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_nam
e "DDX11L1-002"; transcript_source "havana";
1       processed_transcript    exon    11869   12227   .       +       .       gene_id "ENSG00000223972"; transcript_id
 "ENST00000456328"; exon_number "1"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; trans
cript_name "DDX11L1-002"; transcript_source "havana"; exon_id "ENSE00002234944";

And all the genes in the text file have an expression value of 0. What exactly is going wrong in this step? Is the conversion to a SAM file the problem? I'm completely lost so any help would me much appreciated.
Thanks!
EDIT:: For anybody who sees this post cause they have a similar problem. Just note that the newer versions of htseq-count don't require sorted .bam files and, so following the editing of the .gtf file, all I needed to do was convert it to .sam using samtools view -h and then pipe this to htseq-count. Sorting the files prior to this conversion would cause errors.

The obvious first thing to check is that your gtf and genome match. Did you get them from the same place? Specifically, do the chromosome names in the gtf match the ones in the bam?

Good point..... I have just checked, and the gtf file contains no chromosomes actually. After running grep chr over the file it returns nothing. What is the reason for this? And how would I find one that works with the BAM files provided?

"grep chr" won't find anything if the chromosome names are bare numbers.  The question is, does your bam have bare numbers as chromosome names?

@swbarnes2 I converted your answer to a comment since it wasn't giving an answer but only requesting more information from the OP. Please only post answers if they are actually providing an answer. If you need more details, leave a comment instead.

Can you show us a few lines of your gtf file? I am guessing you have 1 instead of chr1.

My .bam file chromosomes are formatted as such: chr1, chr2,....

If your BAM file is using UCSC chromosome names (chr1, chr2, etc.) and your GTF file is using Ensembl chromosome names (1, 2, etc.) then htseq-count won't know how to match your alignments to any of the genes. The simplest solution is to not use htseq-count, but instead to use featureCounts, which will handle the chromosome name difference transparently (it's one of the few programs that does) and is much faster anyway. In the future, download your GTF and fasta files from the same source, don't mix and match them.

Thanks for the reply! My BAM file does indeed use this chromosome name system.Thing is though, my lab wants to continue using htseq for consistency. Would this problem be fixed by editing the .gtf file chromosome names to chr1, chr2 etc. ?

Yes, you can edit the GTF, though if you want real consistency you might see if you were using the exact same fasta and GTF files before.

@h3ab74 you can change the names by simply running sed 's/^/chr/' file.gtf > file.fixed.gtf`` and then use the file.fixed.gtffor your analysis. However, ***and this is important***: make sure that the bam files and the gtf are both referring to the same genome as Devon said, and also make sure there are no blank links in your gtf file. The command will only work correctly if _all_ the lines in the file should start withchr`.

I ran a similar command to that one: sed /^[[:digit:]]/s/^/chr/ file.gtf > newfile.gtf. I checked out newfile.gtf and made sure the chromosome names were corrected; they're now chr1, chr2, etc. However, when I ran the original command with the new file I still got the same problem. Is it maybe an issue with my samtools view command?

Thank you all for your helpful comments, I have made it work. And @DevonRyan, good point, I made sure to use the same .gtf files!


text:
Searching tool to calculate phase/switch error rate
I'm looking for a tool which, given a truth vcf file and a test vcf file, calculates the phase/switch error rate. I performed phasing of a vcf using WhatsHap and want to compare the outcome to some ground truth phased vcf I have. I can't find a tool, and I don't want to write it myself. It seems like a rather common thing... so I'd expect it to exist, but my google-fu is weak today.
Example of test lines:
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample

chr10   100588  .       G       A       .       PASS    AF=0.18;AC=1;NS=150;AN=2;EAS_AF=0.14;EUR_AF=0.39;AFR_AF=0.08;AMR_AF=0.31;SAS_AF=0.2;VT=SNP;DP=11907     GT:PS   1|0:48005
chr10   102385  .       A       G       .       PASS    AF=0.07;AC=1;NS=150;AN=2;EAS_AF=0.02;EUR_AF=0.0;AFR_AF=0.1;AMR_AF=0.19;SAS_AF=0.01;VT=SNP;DP=22665      GT:PS   0|1:48005
chr10   105170  .       G       A       .       PASS    AF=0.03;AC=1;NS=150;AN=2;EAS_AF=0.02;EUR_AF=0.0;AFR_AF=0.0;AMR_AF=0.19;SAS_AF=0.0;VT=SNP;DP=17140       GT:PS   0|1:48005
chr10   105365  .       G       T       .       PASS    AF=0.19;AC=1;NS=150;AN=2;EAS_AF=0.14;EUR_AF=0.39;AFR_AF=0.09;AMR_AF=0.31;SAS_AF=0.2;VT=SNP;DP=21845     GT:PS   1|0:48005
chr10   106057  .       T       C       .       PASS    AF=0.13;AC=1;NS=150;AN=2;EAS_AF=0.02;EUR_AF=0.0;AFR_AF=0.32;AMR_AF=0.19;SAS_AF=0.01;VT=SNP;DP=22559     GT:PS   0|1:48005

Example of truth lines
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample
chr10   100554  .       C       CT,CTT  30      .       .       GT:AD   1|2:0,1,1
chr10   100588  .       G       A       30      .       .       GT:AD   0|1:1,1
chr10   102385  .       A       G       30      .       .       GT:AD   1|0:1,1
chr10   102636  .       T       C       30      .       .       GT:AD   1|1:0,2
chr10   102757  .       T       C       30      .       .       GT:AD   0|1:1,1
chr10   104815  .       GCC     G,GC    30      .       .       GT:AD   2|1:0,1,1
chr10   105170  .       G       A       30      .       .       GT:AD   1|0:1,1
chr10   105365  .       G       T       30      .       .       GT:AD   0|1:1,1
chr10   106057  .       T       C       30      .       .       GT:AD   1|0:1,1
chr10   106110  .       C       G       30      .       .       GT:AD   1|1:0,2
chr10   106261  .       A       G       30      .       .       GT:AD   1|0:1,1
chr10   108612  .       T       G       30      .       .       GT:AD   1|0:1,1
chr10   108646  .       A       AT      30      .       .       GT:AD   0|1:1,1
chr10   108834  .       G       A       30      .       .       GT:AD   1|0:1,1
chr10   110840  .       C       T       30      .       .       GT:AD   1|1:0,2
chr10   111743  .       A       AT      30      .       .       GT:AD   1|0:1,1
chr10   112016  .       T       A       30      .       .       GT:AD   1|0:1,1
chr10   112262  .       C       T       30      .       .       GT:AD   1|1:0,2
chr10   113006  .       C       T       30      .       .       GT:AD   0|1:1,1
chr10   113031  .       G       A       30      .       .       GT:AD   0|1:1,1
chr10   113136  .       G       A       30      .       .       GT:AD   1|0:1,1
chr10   113359  .       A       T       30      .       .       GT:AD   1|1:0,2
chr10   113583  .       C       G       30      .       .       GT:AD   0|1:1,1
chr10   113997  .       G       C       30      .       .       GT:AD   1|1:0,2

Could you give us a few lines of each file to test on? Do you just want to compare the genotype field for the same variant in the two vcf files?

I've edited my question to include example lines. The truth file has every variant phased (it is created based on the synthetic diploid method https://www.nature.com/articles/s41592-018-0054-7 and https://github.com/lh3/CHM-eval/tree/master/dip-call).

Assuming your files are as you show, and the genotype field is always the 1st entry of the sample information field (the 10th), then you can just do this directly with a simple awk script:
$ awk -F"\t" -vOFS="\t" 'BEGIN{print "#CHROM","POS","ID","REF","ALT","Truth","Test"}{ if(/^[^#]/){var=$1$2$4$5; gt=substr($10,1,3); if(NR==FNR){a[var]=gt; next;} if(var in a && a[var]!=gt){print $1,$2,$3,$4,$5,gt,a[var]}}}' truth.vcf test.vcf
#CHROM  POS     ID  REF ALT Truth   Test
chr10   100588  .   G   A   1|0 0|1
chr10   102385  .   A   G   0|1 1|0
chr10   105170  .   G   A   0|1 1|0
chr10   105365  .   G   T   1|0 0|1
chr10   106057  .   T   C   0|1 1|0

Here's the script split into separate lines for clarity (you can still copy paste directly into a terminal):
awk -F"\t" -vOFS="\t" '
    BEGIN{
        print "#CHROM","POS","ID","REF","ALT","Truth","Test"
    }
    {
        if(/^[^#]/){
            var=$1$2$4$5;
            gt=substr($10,1,3);
            if(NR==FNR){
                a[var]=gt;
                next;
            }
            if(var in a && a[var]!=gt){
                print $1,$2,$3,$4,$5,gt,a[var]
            }
        }
    }' truth.vcf test.vcf

The basic idea is:

awk -F"\t" -vOFS="\t": set the input and output field separators to a tab.
BEGIN{print "#CHROM ... }: print a header.
if(/^[^#]/){: skip the header lines of the input vcf files
var=$1$2$4$5: get a "signature" for this variant, by concatenating the chromosome, position, ref and alt fields.
gt=substr($10,1,3); : get the genotype value.
if(NR==FNR){: if this is the 1st file/
a[var]=gt; next: save the genotype for this variant and skip to the next line.
if(var in a && a[var]!=gt){: this part will only be run for the 2nd file because of the 'next' in the previous block. If this variant was in the truth set and if its genotype there isn't the same as in the test set, print the variant and the two genotypes.

If you need the error rate, you can extend the script thusly:
awk -F"\t" -vOFS="\t" '
BEGIN{
    print "#CHROM","POS","ID","REF","ALT","Truth","Test"
}
{
    if(/^[^#]/){
        var=$1$2$4$5;
        gt=substr($10,1,3);
        if(NR==FNR){
            a[var]=gt;
            next;
        }
        tot++;
        if(var in a && a[var]!=gt){
            bad++;
            print $1,$2,$3,$4,$5,gt,a[var]
        }
    }
}
END{
    printf "Correct: %s of %s (%s%)\n",tot-bad,tot,(tot-bad)*100/tot;
    printf "Wrong: %s of %s (%s%)\n",bad,tot,bad*100/tot;
}' truth.vcf test.vcf

On your example dataset, that prints:
#CHROM  POS ID  REF ALT Truth   Test
chr10   100588  .   G   A   1|0 0|1
chr10   102385  .   A   G   0|1 1|0
chr10   105170  .   G   A   0|1 1|0
chr10   105365  .   G   T   1|0 0|1
chr10   106057  .   T   C   0|1 1|0
Correct: 0 of 5 (0%)
Wrong: 5 of 5 (100%)

In the test VCF, "PS" means "phase set". Switches between phase sets are not counted as errors. In general, we do need a proper script for evaluation. There is not a common script probably because the inputs often change a lot.

@user172818 oh wow, thanks! I hadn't encountered the PS before. Does this mean my simplistic approach is basically useless?

I think it will be useful when "PS" is not involved. For example, Hi-C phasing is supposed to give whole-chromosome phasing. We don't need "PS". Your script is applicable in that case.

Yes, while the truth set has "left" and "right" haplotypes (paternal or maternal) the test set is split in haplotype/phase blocks. Within a block, we can figure out which SNPs are on the same haplotype (or aren't), but not across blocks.  My example is not perfect, but I believe it's 100% correctly phased here.


text:
Differentially methylated position analysis in a related sample?
I'm trying to figure out how to do a DMP analysis (using minfi dmpFinder) on a related sample (if it's even possible). Right now the code (not written by me) is:
dmps = dmpFinder(mVals,
                 pheno = targets1$X8yrfactor,
                 type = "continuous")
pheno.exists <- !is.na(targets1$X8yrfactor)
pheno <- targets1$X8yrfactor[pheno.exists]
mVal.phenoNotNA <- mVals[, pheno.exists]
dmps = dmpFinder(mVal.phenoNotNA, pheno = pheno, type = "continuous")
dmps$probe = paste(rownames(dmps))
rownames(dmps) = NULL
dmps = dmps[,c(6,1:5)]
write.csv(dmps, file = "dmps_8yr_minfi.csv", row.names = F)

The sample size is 96, i.e. 48 twin pairs, which is why we want to be able to control for relatedness in the sample. By "twin pairs", I mean the methylation data we're using is from a sample of 48 pairs of identical human twins. At the moment I believe we're only going to be testing one at a time if we do find a way to test them on this sample.
Is there anyway to code it similarly to a mixed linear model where a randomization variable for each related individual is taken into account? I've been referred/found this, this and this resource but I'm not sure they capture what I need. Correct me if I'm wrong.
A simple linear model will not do this, which is why we're trying to figure out how to do this similarly to a mixed model. 'dmpFinder' identifies CpGs where methylation is associated with a continuous or categorical phenotype ('pheno', which we have multiple we want to test, we're trying to write the code first). 'targets1' is our cases (3 of them) that were removed following QC.

If you have multiple phenotype you need to control by that and consider the effect of one category into other categories or interaction between categories. Could you post the header of the phenotype (all the phenotype you want to test and what it is)?


text:
extracting gene names from pubmed abstracts
I have been extracted number of abstracts from pubmed. now I am having a bunch of abstracts of pubmed articles now from this abstracts I need to extract gene names which are included in abstracts.

What you need is a list of genes a priori. In my opinion the problem is there are often no "markers" to alert a text finding algorithm that XXX is gene or a protein, i.e. often no labelled as gene or protein. I could be wrong however.

Welcome to Bioinformatics.StackExchange @VenkataraoMannem! Please take a look at our Help Center, in particular how to ask a good question. Specifically, please tell us what you have tried so far and why it is not working. There is a lot of published literature on text mining of biomedical literature. Are you having trouble finding this literature? (Have you even made an attempt?) Have you found a strategy but now need help running it?  Please use the [edit] button to update your post with these details so we can better help you!

Cross-posted from Biostars: https://www.biostars.org/p/368251/#368431


text:
Variant calling without matched normal sample
I have WGS .bam files for 3 patients (tumour and its matched derived model namely organoid) but I don't matched normal sample. If I call variants of each patients (tumour and its matched organoid), how I can use read counts at germline heterozygous positions estimated by GATK 3.2-2 to compensate for the absence of matched normal sample? I heard people use dbsnp VCF instead of the matched normal small variant VCF. But, I don't know start from where? which GATK function does that?

I should mention, Calling SNV and indel in many tools returns vcf but
  I called copy number by varscan that did not return .vcf output so I
  am not sure what to do for CNV

Any helps please?

If you just want to filter out calls present in dbSNP then use:
java -jar GenomeAnalysisTK.jar \
   -T SelectVariants \
   -R reference.fasta \
   -V patient.vcf \
   --discordance dbSNP.vcf \
   -o patient.filtered.vcf

--discordance will produce calls not present in dbSNP.

Thanks a lot @Devon, Calling SNV and indel in many tools returns vcf but I called copy number by varscan that did not return .vcf output, is this command is applicable yet?

It's usually best to mention things like that in your post... I presume this will still work, but I can't say I've ever tried.

https://software.broadinstitute.org/gatk/documentation/article?id=11682

I heard GenomeAnalysisTK.jar has been replaced. Do you know by which function we can filter a vcf file based on dbSNP?

Here is a clue by scatngs
https://github.com/cancerit/ascatNgs/wiki/Human-reference-files-from-1000-genomes-VCFs
and already calculated file is here
ftp://ftp.sanger.ac.uk/pub/cancer/support-files/CPIB/ascatNgs/Human/GRCh37/
Although installing scatngs itself is a big burden


text:
Finding genes especific to microenvironment
I have RNA-seq .bam files for 3 patients, tumour and its matched derived model, namely organoid, but I don't have any matched normal sample.
Differentially expressed genes between a tumour and its model mainly come from the presence and lack of microenvironment in cancers and organoids respectively and response to stimuli (artificial environment in organoid culture vs cancer.
How I can get the differentially expressed genes specific to microenvironment?
Any intuition please?

Sauret package for unsupervised learning is a start, lots of discussion on this board about it

I can't find a package named Sauret, if you refer to Seurat it is for scRNA-seq which this question is not tagged with it so I think it won't work well in this case.

The genes related to the environment can be identified using a linear model where each organoid is compared to the correspondent original tumor. This is assuming that the environment for all the organoids is the same and that the differences between tumors is just from different patient but are from the same organ/region.


text:
I am getting a parsing error when uploading .chip file to GSEA software
I get the following error when I try to load a .chip file to GSEA software. What is the problem?
---- Full Error Message ----
There were errors: ERROR(S) #:1
Parsing trouble
java.lang.NumberFormatException: ...

---- Stack Trace ----
# of exceptions: 1
------For input string: "YBL028C"------
java.lang.NumberFormatException: For input string: "YBL028C"
    at sun.misc.FloatingDecimal.readJavaFormatString(Unknown Source)
    at sun.misc.FloatingDecimal.parseFloat(Unknown Source)
    at java.lang.Float.parseFloat(Unknown Source)
    at edu.mit.broad.genome.parsers.TxtDatasetParser._parseNoDesc(TxtDatasetParser.java:169)
    at edu.mit.broad.genome.parsers.TxtDatasetParser.parse(TxtDatasetParser.java:131)
    at edu.mit.broad.genome.parsers.TxtDatasetParser.parse(TxtDatasetParser.java:87)
    at edu.mit.broad.genome.parsers.ParserFactory.readDatasetTXT(ParserFactory.java:202)
    at edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:749)
    at edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:726)
    at edu.mit.broad.genome.parsers.ParserWorker.doInBackground(ParserWorker.java:52)
    at javax.swing.SwingWorker$1.call(Unknown Source)
    at java.util.concurrent.FutureTask.run(Unknown Source)
    at javax.swing.SwingWorker.run(Unknown Source)
    at java.util.concurrent.ThreadPoolExecutor.runWorker(Unknown Source)
    at java.util.concurrent.ThreadPoolExecutor$Worker.run(Unknown Source)
    at java.lang.Thread.run(Unknown Source)

Edit:

I am also using GSEA 3.0. What does not make sense is the other files have the right extensions and are not getting this error, but only this file is failing to upload.

https://www.gsea-msigdb.org/gsea/index.jsp
"Gene Set Enrichment Analysis (GSEA) is a computational method that determines whether an a priori defined set of genes shows statistically
significant, concordant differences between two biological states
(e.g. phenotypes). "

Hi there, would you mind sharing a bit more details? What have you tried to resolve your problem? DO you think it's a problem of installation or the file? Note also that you can [edit] your question.

Also perhaps it could be useful to know which version of GSEA are using and how are the first lines from the file (or those around YBL028C )

Please show us an example of the file you are trying to load. It is complaining about YBL028C. That looks like a gene/protein name, but the tool is trying to read it as a number. So where do you have YBL028C in your data? What is a .chip file?

did you already solve the problem? I also came across the same error when uploading a similar .chip file.