The problems in the design of data processing pipelines

The design problem is the most big question for me, because this will determine whether your work will be useful in next a few days. It's much important then whether you use Emacs or some automatic unit test harness. After struggling these days, I find some pitfalls in my design process.

 

Design when programming

For a huge pipeline, I tend to design the next step when programming this step. This frequently leads to doing a lot of unnecessary work. On the other hand, we can't design all the steps before start. But this should not be the excuse of the avoid design early. [Solution], design as much as possible.

Optimize too early
Last week, I took two days to learn a utility in GATK that can accelerate the speed of retrieving fasta sequence of certain interval from whole genome.  I had done this in R, but it's very slow. Now, when the pipeline is almost done, I find that I don't use accelerate utility at all, and the speed of R is acceptable. So the two days I spent on this utility is kind of wasting time.So the solution of this design trap: Do the necessary things, not the beautiful things.
 

Paralyzed in front of multiple choices

When facing a lot of choices, my brain is paralyzed. And then I tried to escape from such situation, and do some distracting things. In such situation, write down sub steps of each choice. Then you'll get a feeling of which one you like, but not which one is better.

vcf-compare problems in practice

Recently, I'm using vcftool to compare two vcf files. The command to finish this job is vcf-compare. But the developer of vcf tools doesn't illustrate the how to use it in detail, which causes a lot of problem in pratice.

The two vcf files I'm comparing are supposed to overlap a lot with each other. But the vcf compare told me that there was no overlap between them. This is the output of vcf-compare:

Error: There is no overlap between any of the samples, yet haplotype comparison was requested.

By chance, there should be some overlapping. So I supposed there were some problems in my pipeline and data. After comparing a lot of situations, I found two problems of the data:
1. The format header should be same. Since the two vcf files come from different bam files, the last field of format header is the name of their original bam file. This will lead to the difference in format and it should be same.

2. The name of chromosome. In my case, one file using 1 to representing chr1, while the other using chr1. So I need to add chr in before the numbers.

After fixing two problems above, vcf-compare works well.

Emacs Mark Region and Scrolling Problem

The problem looks like this in emacs when I set the mark:

1. C+SPC set the mark at the beginning of the paragraph

2. Move to the end to paragraph.

3. If I need to scroll the page in step 2, the cursor will move to the current view of the emacs

4. Then I use left click to select the region between the start mark and current cursor. But I only get part of the region that in the current view.

After making me depressed a long time, I realized that the reason is due to the failing of set mark in step 1. C+SPC is not captured by Emacs, but by my Input Method.

Synchronize the folder to the server in the Rstudio

1. Install the WinSCP

2. Set WinSCP installation path to the WIN7 Environment Variable in $PATH

3. Write the code:
1. Winscp_sync.txt:

                  open stored_session_name_in_WinSCP

                  synchronize both E:\Projects\R\data\server /homed/home/shi/anthony/tfbs_pwm/rsnp

                  exit

2. Winscp_sync.bat:

                WinSCP.exe /console /script=winscp_sync.txt

4. In Rstudio

               Call: system('./winscp_sync.bat')#Path needed

Matlab with large input data

Yesterday, I was dealing with a big file (about 600MB). At first, I tried to use "textread" to input it to matlab. Though my computer equiped with 8GB physical memory, the matlab finally ended in a long list of error.

Then, I have to use the a much more slower method, read line by line. "fgetl" and "sread", the read the string and the second parse the string according to the given format, just like scanf in C. However, the final program run very slow, about 1440s to finish.

At last, I fould another function 'textscan', which is very similar to textread, but more powerful. It can be use

C = textscan(fid, 'format', N) 

It can read the data part by part of a scale in N, which strikes a better balance between the textread and fgetl. In this way, the operation is paralleled to some extent, so the speed would be faster.

By the way, the awk is a poweful tool to deal with text.

Manage the Complexity

The essential of software is to manage the complexity, because the complexity of one idea to implantation  grow with a factor of 10^9. If there is no rule, the second thermal law will work silently. And your work will be in a huge mess.

And we can also see the complexity in research, daily life, and ...