last day (18 days later) » 

10:59 AM
1
A: Count the number of times a substring appears in a file and place it in a new column

Ted LyngmoSince grep -c seems to give you the correct count (matching lines, not counting multiple occurances on the same line) you could read the 7 fields from the TSV (BED) file and just print them again with the grep output added to the end: #!/bin/bash # read the fields into the array `v`: while read ...

 
Hi Ted, just gave this a go and works! Seen you've just edited it so I'll give that a run, it does exactly what I want, but I think the size of the files I'm working with means it might be running a touch slow!
 
@KLM117 Great! Regarding speed: Yeah, it'll be running grep (opening and closing the fasta file) one time per line in the tsv.bed file so speed could become an issue. You could perhaps use readarray to read the whole fasta file into memory one time and then find your entries in that array. Might be faster or it may not be. Not sure.
@KLM117 Just for fun, I rewrote this in C++ using a memory mapped file on "hg19.fasta" and like to compare it with the above bash script. How many lines are there in tsv.bed and hg19.fasta? If they don't contain anything copyrighted I wouldn't mind downloading the actual files.
 
that awesome! So the hg19.fasta is the human reference genome 19, this is relatively massive file with 61,998,735 lines (~3gb). This is publically available and I tend to ftp it from the UCSC hgdownload.soe.ucsc.edu/downloads.html using this code: ftp hgdownload.soe.ucsc.edu #use "Anonymous" as name, and email as passwd. cd goldenPath/hg19/bigZips mget hg19.fa.gz gzip hg19.fa the tsv file tends to be around 600,000 lines (i'm actually running it on around 12 tsv files) I generated them myself so happy to send you one to try out with!
 
@KLM117 I got hg19.fa. Thanks! 4796636 lines contain all N:s. Perhaps one could treat them specially? Would it be ok to just ignore those lines or will you do matching against those too? And yes please if you could put one or tsv files up somewhere that'd be grand! One strange thing. Your first pattern, ATGGTGCTTTGTTATGGCAGCTC isn't found in hg19.fa anywhere ...
 
Hi Ted, yes, I think that removing the N's could be removed, they simply stand for any nucleotide (A,C,G, or T) and so aren't needed. I've thrown an example tsv in a github repository github.com/kai-lawsonmcdowall/tsv_file, it's ~420,000 lines/25mb. Let me know how it goes!
 
10:59 AM
@KLM117 Nice! What's bothering me is that many of the entries in the tsv file are not found in the geonome file and I've found no entries with more than 1 match in ~60 first entries. Is that how it should be?
 
Hey Ted, I think this comment is moved to chat to prevent a super-long string of comments, hope that is ok! Yes, or at the very least it's ok if there aren't any matches in the genome, these are predicted sequences so they simply may not exist in the hg19.
 
 
3 hours later…
1:41 PM
@KLM117 Aha, ok. Ok, I have just been away and had my third covid-shot and started looking it this again. It's looking really good. I've made it scale so that all CPU cores run on max when searching. I'm running on an old server with only 12 threads and a spinning disk. With an SSD disk and more threads it'd probably be a lot faster. I'll try it on my Windows-machine later to see how that goes.
@KLM117 You don't happen to have development tools for C++ already installed? Like a C++20 compiler + the boost library?
 
2:19 PM
@TedLyngmo, Ah very good, I got mine on monday (hope you are symptom-free)! I'm very glad to hear that. Currently I'm running 8 cpus/ssd I've got access to AWS ec2 instances, so I could in theory run a 48 vCPU/ssd instance. I believe I do. Again thank you, I really appreciate this!
 
Cool - that should speed things up a lot . I'll continue tinkering with it. I don't think it'll be appropriate for an answer to the question at Stackoverflow but I'll put it on github when I have something ready and we can take it from there.
 
2:50 PM
Awesome, keep me in the loop and let me know how it goes!
 
 
1 hour later…
3:53 PM
@KLM117 Ok, I have a first draft up at github.com/TedLyngmo/hgsearch - I haven't tested it on my fast machine with lots of cores + SSD but it's a whole lot faster than my bash script even on my old server :-)
 
 
3 hours later…
7:07 PM
@KLM117 Some initial tests on my faster machine (not with Visual Studio but with g++ in WSL) shows that it's a lot faster there but I can't help to wonder if it wouldn't be even faster with VS. Curretly it seems to process 4-5 of your tsv entries per second against hg19.fa (with the NNNNNN-lines removed). That means it'll take more than a day to process the whole file.
@KLM117 I also see that upper and lowercase being mixed in hg19.fa. What does that mean?
 
 
3 hours later…
10:03 PM
@TedLyngmo that's great, I don't have the files right now but I shall take a look first thing tomorrow morning! In terms of the upper and lower case letters, lower case letters are something called soft-masked reads. They are essentially useful when parts of the genome (the fasta file) are highly repetitive or low complexity. does the C++ also search for lowercase or exact match? If possible I think only exact matches makes more sense
@TedLyngmo Re the g++/wsl, I think the amazon Linux essentially reflects this. Would running with VS mean running the code from it's API? Out of curiosity, how powerful is your faster machine?
 
10:49 PM
@KLM117 "upper and lower case" - Ok, great, the program does exact matches right now. "g++/wsl" - Aha, ok, I'd expect them to run linux container on a linux host - but I don't think it'll matter that much, It'll be fine.
It's just that I've noticed a brutal performance when running directly on the WinAPI with the execution policies I put in similar programs in the past, but the g++ performance isn't bad at all.
I hope I can test it later. Right now I've only made some comparisons between the different solutions to your problem (stored in stats/bed in the repo). bed 1 means just running 1 of your lines in the tsv-bed file, 2, two lines etc.
My faster machine is a Core i9 with 12 cores, 24 hyperthreads. It doesn't seem to spin up all CPU:s to max frequency even though the program runs all 24 threads on 100% CPU :-) I hope that's something in WSL. it'd be really nice to get it do run on max CPU frequency. :)
 

  last day (18 days later) »