# Missing records in my FASTA file when writing new FASTAs

**URL:** <https://discourse.julialang.org/t/missing-records-in-my-fasta-file-when-writing-new-fastas/74032>\
**Category:** Biology, Health, and Medicine\
**Tags:** question\
**Created:** [January 4, 2022, 2:15pm UTC](https://discourse.julialang.org/t/missing-records-in-my-fasta-file-when-writing-new-fastas/74032 "2022-01-04T14:15:35Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Lamma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lamma/32/31006_2.png) [@Lamma](https://discourse.julialang.org/u/Lamma)\
**Post date:** [January 4, 2022, 2:15pm UTC](https://discourse.julialang.org/t/missing-records-in-my-fasta-file-when-writing-new-fastas/74032/1 "2022-01-04T14:15:35Z")

</div>

I have a co-assembly and some alignment data telling me if a reads from a sample that helped spawn that co-assembly align to a contig. I am using this to write a fasta file for that sample.

```julia
    t1 = Threads.@spawn DataFrame(CSV.File(depth, delim="\t"))
    abundf = fetch(t1)
    rename!(abundf, [:Contig,:trimmed_mean])

    present = @subset!(abundf, :trimmed_mean .!= 0)

    passContig = Set(present[:,"Contig"])

```

Here I am reading in and removing all rows that have a depth of 0 indicating that no reads from that sample aligned to that contig (where each row is a contig in the co assembly)

```julia
reader = FASTA.Reader(GzipDecompressorStream(open(assembly))) #if gzipped
    # Generating the writer to push the output into
    writer = open(FASTA.Writer, "D:/OneDrive - University of Copenhagen/PhD/Projects/Termite_metagenomes/assembly/matriline/bySample/"*sample*".fasta")
    println("Writing to file")
    # Looping over the contigs
    record = FASTA.Record()
    while !eof(reader)
        read!(reader, record)
        if in(FASTA.identifier(record), passContig)
            write(writer, record)
        end
    end

    close(reader)

```

Next I am opening the co-assembly and going over its headers and seeing if that header is in the `passContig` list for the given sample.

However the length of `passContig` is not equal to the resulting number of reads in the fasta being written. This should not be the case as any contig present in the `passContig` is in the co-assembly.

```julia
> length(passContig)
121964

```

whilst the number of reads in the fasta is: 121955 (generated from `grep -c ">" sample.fasta`)

Therefore I next tried to identify what contigs were being missed:

```julia
reader = FASTA.Reader(GzipDecompressorStream(open(assembly))) #if gzipped
headers= Set([FASTA.identifier(rec) for rec in reader])
close(reader)

for r in passContig
    if !(r in headers)
        println(r)
    elseif r in headers
        println("it is here")
    end
end

```

However this results in nothing but “it is here”. Would anyone know why this is happening?

---

<div class="post-metadata">

**Author:** ![kevbonham](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kevbonham/32/216165_2.png) [@kevbonham](https://discourse.julialang.org/u/kevbonham)\
**Post date:** [January 4, 2022, 4:49pm UTC](https://discourse.julialang.org/t/missing-records-in-my-fasta-file-when-writing-new-fastas/74032/2 "2022-01-04T16:49:59Z")

</div>

Not knowing how the CSV file was generated, my only guess is that it has some headers that aren’t in your fasta. Was there a trimming or filtering step that happened after the CSV was generated, but before you’re doing this filtering?

If it were me debugging, I’d do a `setdiff()` to find which headers are in the dataframe but not in the fasta file, then start digging around in the file.

---

<div class="post-metadata">

**Author:** ![Lamma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lamma/32/31006_2.png) [@Lamma](https://discourse.julialang.org/u/Lamma)\
**Post date:** [January 5, 2022, 8:41am UTC](https://discourse.julialang.org/t/missing-records-in-my-fasta-file-when-writing-new-fastas/74032/3 "2022-01-05T08:41:52Z")

</div>

I did many variation of checking (including how you suggested) but in the end it was because I did close the writer!!

---

<div class="post-metadata">

**Author:** ![jonathanBieler](https://avatars.discourse-cdn.com/v4/letter/j/82dd89/32.png) [@jonathanBieler](https://discourse.julialang.org/u/jonathanBieler)\
**Post date:** [January 5, 2022, 10:27am UTC](https://discourse.julialang.org/t/missing-records-in-my-fasta-file-when-writing-new-fastas/74032/4 "2022-01-05T10:27:36Z")

</div>

As a side note you could give [BioRecordsProcessing](https://github.com/jonathanBieler/BioRecordsProcessing.jl) a try, I wrote it exactly to reduce the boilerplate when doing this kind of file operations :

```julia
out_path = "D:/OneDrive - University of Copenhagen/PhD/Projects/Termite_metagenomes/assembly/matriline/bySample/"

BioRecordsProcessing.process(FASTX.FASTA, assembly, out_path) do record
    in(FASTA.identifier(record), passContig) && return record
    return nothing
end

```