CGI/Perl Guide | Learning Center | Forums | Advertise | Login
Site Search: in

  Main Index MAIN
INDEX
Search Posts SEARCH
POSTS
Who's Online WHO'S
ONLINE
Log in LOG
IN

Home: Perl Programming Help: Intermediate: Re: [Chris Charley] DNA Sequence Count Perl Program HELP!: Edit Log



Laurent_R
Veteran / Moderator

Oct 28, 2012, 3:02 AM


Views: 6406
Re: [Chris Charley] DNA Sequence Count Perl Program HELP!

Hi,

I just made a test with Chris Charley's DNA sequence file. Contrary to what I thought and wrote earlier, the end of line characters alter the computed skew quite significantly.

Example of a few windows on the DNA file as posted by Chris:


Code
Pass 1 Window Number 1 : there are 424 guanine and 443 cytosine in this window. Skew value is: -0.0989374637120816 
Pass 1 Window Number 2 : there are 451 guanine and 369 cytosine in this window. Skew value is: 0.542005420054201
Pass 1 Window Number 3 : there are 363 guanine and 359 cytosine in this window. Skew value is: 0.0308644356823741


I then removed the new line/carriage return characters from the same file:


Code
$ perl -i -pe 's/\r\n//g' test.txt


And, running the same test with the new file without the carriage returns:


Code
Pass 1 Window Number 1 : there are 437 guanine and 469 cytosine in this window. Skew value is: -0.150618713433777 
Pass 1 Window Number 2 : there are 467 guanine and 375 cytosine in this window. Skew value is: 0.582739509105305
Pass 1 Window Number 3 : there are 395 guanine and 389 cytosine in this window. Skew value is: 0.0393473584806673


So, Kevone, you should remove these new line/carriage return characters from your file if you are going to use my suggested program (you can do it with my Perl one-liner above, but you may have to change it depending on whether you are on Unix or Windows).

Another point is that I found a bug in my calculation of the skew: in the formula, I use $window_size for the total number of nucleobases analyzed, whereas I should have used the actual number of bytes read, otherwise the calculated skew is wrong whenever I reach the end of the file (i.e. for each last window of each pass).

So the line calculating the skew is now as follows:


Code
		my $skew_value = ($size_read / $nr_C)*($nr_G - $nr_C)/($nr_G + $nr_C);


Without that change, I had the following result for the first pass:



Code
Pass 1 Window Number 1 : there are 437 guanine and 469 cytosine in this window. Skew value is: -0.150618713433777 
Pass 1 Window Number 2 : there are 467 guanine and 375 cytosine in this window. Skew value is: 0.582739509105305
Pass 1 Window Number 3 : there are 395 guanine and 389 cytosine in this window. Skew value is: 0.0393473584806673
Pass 1 Window Number 4 : there are 467 guanine and 416 cytosine in this window. Skew value is: 0.277680982663995
Pass 1 Window Number 5 : there are 486 guanine and 436 cytosine in this window. Skew value is: 0.248761169376505
Pass 1 Window Number 6 : there are 509 guanine and 506 cytosine in this window. Skew value is: 0.0116824704530851
Pass 1 Window Number 7 : there are 546 guanine and 534 cytosine in this window. Skew value is: 0.0416146483562214
Pass 1 Window Number 8 : there are 552 guanine and 524 cytosine in this window. Skew value is: 0.0993217741706632
Pass 1 Window Number 9 : there are 507 guanine and 555 cytosine in this window. Skew value is: -0.162874739145926
Pass 1 Window Number 10 : there are 194 guanine and 217 cytosine in this window. Skew value is: -0.515770235572449
Pass 2 Window Number 1 : there are 421 guanine and 482 cytosine in this window. Skew value is: -0.280301254922504
Pass 2 Window Number 2 : there are 454 guanine and 358 cytosine in this window. Skew value is: 0.660483804386713


With this correction, I have now this:


Code
Pass 1 Window Number 1 : there are 437 guanine and 469 cytosine in this window. Skew value is: -0.150618713433777 
Pass 1 Window Number 2 : there are 467 guanine and 375 cytosine in this window. Skew value is: 0.582739509105305
Pass 1 Window Number 3 : there are 395 guanine and 389 cytosine in this window. Skew value is: 0.0393473584806673
Pass 1 Window Number 4 : there are 467 guanine and 416 cytosine in this window. Skew value is: 0.277680982663995
Pass 1 Window Number 5 : there are 486 guanine and 436 cytosine in this window. Skew value is: 0.248761169376505
Pass 1 Window Number 6 : there are 509 guanine and 506 cytosine in this window. Skew value is: 0.0116824704530851
Pass 1 Window Number 7 : there are 546 guanine and 534 cytosine in this window. Skew value is: 0.0416146483562214
Pass 1 Window Number 8 : there are 552 guanine and 524 cytosine in this window. Skew value is: 0.0993217741706632
Pass 1 Window Number 9 : there are 507 guanine and 555 cytosine in this window. Skew value is: -0.162874739145926
Pass 1 Window Number 10 : there are 194 guanine and 217 cytosine in this window. Skew value is: -0.219202350118291
Pass 2 Window Number 1 : there are 421 guanine and 482 cytosine in this window. Skew value is: -0.280301254922504
Pass 2 Window Number 2 : there are 454 guanine and 358 cytosine in this window. Skew value is: 0.660483804386713


As it can be seen, the computed skew values are the same for every window, except the last window of each pass (window number 10 in the above example). So I think you should probably apply the same correction.

EDIT: I checked Chris Charley's program, it has the same bug in the skew calculation (since Chris probably reused my formula). So you should modify his program if it is the one you are going to use.


(This post was edited by Laurent_R on Oct 28, 2012, 3:08 AM)


Edit Log:
Post edited by Laurent_R (Veteran) on Oct 28, 2012, 3:08 AM


Search for (options) Powered by Gossamer Forum v.1.2.0

Web Applications & Managed Hosting Powered by Gossamer Threads
Visit our Mailing List Archives