Originally posted by reubennowell
View Post
Seqanswers Leaderboard Ad
Collapse
Announcement
Collapse
No announcement yet.
X
-
I found another case where this was happening I thought I'd just add to the list for future reference.
My input fasta file was line wrapped, and there was a blank line between sequences. The problem sequence was where coincidentally the seq ended exactly at the line wrap position, ie:
Code:>not_a_problem_seq AATCGACGTACGTAGCTGATC ATCGATGCTAGCTATAATCGT ACTAGCTACTGG >problem_seq ACTGCTAGCTGATCGATCGTAG ACGTACGTAGCTAGCTGACTGA ACTGATCGTAGCTAGCTGATCG <- here! >next_seq ATCGACGTAGCTAGCTAGCTAT ACTGATGCTACTAGCTGATGCT ACTGA
Leave a comment:
-
In my case, the actual problem was different (I found out later).
I extracted the genomic sequences from UCSC table browser in fasta format and saved it as my reference file. Invoking samtools faidx was giving the error 'different line length...'. However, fixing empty lines did not stop the error, and I discovered that my fasta file was truncated due to timeout while downloading from UCSC; the error message was the last lines of my fasta file and it was causing samtools emit this message.
Leave a comment:
-
That works too. But there are canned answers to this sort of question nowadays, fortunately, that didn't exist in 2009 when I originally posted. Picard, for example, has a FASTA-fixer:
I suspect BioPerl and similar packages were doing this even in 2009, but I was pretty new to the field then and didn't know about them.
Leave a comment:
-
Thanks.
I did this:Code:awk '(NF!=0){print}' sample.fa>tmp.fa awk '(NF==0){print "\blank"}' tmp.fa # check mv tmp.fa sample.fa
Leave a comment:
-
Turns out that in my case it's actually a blank line that's causing the problem. I started wondering why the issue was coming up at all, given that awk says the longest line in the file is 70 characters long:
awk '{ if (x < length()) x = length() } END { print x }' hsap_36.1_hg18.fa
I thought surely it couldn't be complaining about short lines...so I looked for sub-70-character lines:
awk 'length($0) < 70' hsap_36.1_hg18.fa
I got the expected set of header lines, and one trailing incomplete line per header block. But the last line of the chr13 block was blank due to the actual final line of sequence being exactly 70 characters.
So the brute-force solution for me is just to eliminate blank lines from the file:
grep . hsap_36.1_hg18.fa > tmp; mv tmp hsap_36.1_hg18.fa
This still leaves the questions of (1) why samtools can't handle a simple silly blank line and (2) how to fold a FASTA file that does have anomalous line lengths...but my immediate crisis is over. Thanks in advance, though, for any insights anyone has on either of those questions!
Leave a comment:
-
Samtools fasta line length segfault, redux
There's been at least one thread on this already, so I'm wondering if anyone has a handy shell script they could share for fixing .fa files with bad line lengths. In my particular case, I get this error from samtools faidx:
[fai_build_core] different line length in sequence 'chr13'.
...when using our standard reference file, hsap_36.1_hg18.fa.
I want to reformat this file rather than downloading some known-correct one, since I don't know exactly where our file came from and what differences there might be between it and another.
The thread linked above has a snippet of BioPerl, but alas, I don't have that available on the server where I'm working. The last poster outlines an approach, but doesn't provide actual code.
Hopefully I'm not being too lazy here by asking if someone's got such a script in their back pocket...I'm just a bit overwhelmed at the moment. I would love not to have to solve this particular problem from scratch, and I'm sure a solution posted here would be useful to others in the future. Thanks!Tags: None
Latest Articles
Collapse
-
by seqadmin
The complexity of cancer is clearly demonstrated in the diverse ecosystem of the tumor microenvironment (TME). The TME is made up of numerous cell types and its development begins with the changes that happen during oncogenesis. “Genomic mutations, copy number changes, epigenetic alterations, and alternative gene expression occur to varying degrees within the affected tumor cells,” explained Andrea O’Hara, Ph.D., Strategic Technical Specialist at Azenta. “As...-
Channel: Articles
07-08-2024, 03:19 PM -
ad_right_rmr
Collapse
News
Collapse
Topics | Statistics | Last Post | ||
---|---|---|---|---|
Started by seqadmin, Yesterday, 06:46 AM
|
0 responses
9 views
0 likes
|
Last Post
by seqadmin
Yesterday, 06:46 AM
|
||
Started by seqadmin, 07-24-2024, 11:09 AM
|
0 responses
26 views
0 likes
|
Last Post
by seqadmin
07-24-2024, 11:09 AM
|
||
Started by seqadmin, 07-19-2024, 07:20 AM
|
0 responses
160 views
0 likes
|
Last Post
by seqadmin
07-19-2024, 07:20 AM
|
||
Started by seqadmin, 07-16-2024, 05:49 AM
|
0 responses
127 views
0 likes
|
Last Post
by seqadmin
07-16-2024, 05:49 AM
|
Leave a comment: