Counting reads in Cv gonad data

[yaaminiv]

I want to use a loop to count reads in my raw C. virginica gonad sequence data. I tried a loop in this Jupyter notebook but got a gunzip error:

I've tried the echo $(( $(gunzip -c FILENAME | wc -l) / 4 )) code on a single file so I know it works...just not sure how the loop changes things.

[sr320]

Why do you want to count reads? This is provided in Bismarck thanks Steven Steven B. Roberts Kenneth K. Chew Endowed Professor School of Aquatic and Fishery Sciences University of Washington, Seattle WA sr320@uw.edu - 206.866.5141 robertslab.info - @sr320

…

On Mar 18, 2019, 8:54 AM -0700, Yaamini Venkataraman ***@***.***>, wrote: I want to use a loop to count reads in my raw C. virginica gonad sequence data. I tried a loop in this Jupyter notebook but got a gunzip error: I've tried the echo $(( $(gunzip -c FILENAME | wc -l) / 4 )) code on a single file so I know it works...just not sure how the loop changes things. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or mute the thread.

[yaaminiv]

Bismark has trimmed reads, not raw reads. I thought it would be good to include read counts from both untrimmed and trimmed files.

[kubu4]

You're missing the * needed as a wildcard to identify which files you want your for loop to operate on.

[sr320]

True - FastQC should have those numbers

…

On Mar 18, 2019, 9:03 AM -0700, Yaamini Venkataraman ***@***.***>, wrote: Bismark has trimmed reads, not raw reads. I thought it would be good to include read counts from both untrimmed and trimmed files. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[yaaminiv]

@kubu4 RIP me. Thanks!

@sr320 Got it. Total sequences = total reads

[yaaminiv]

Another loop question...

I was able to successfully count reads, but I want to save the grep output in a new file. When I try redirecting the output to a new file, I only get the last entry saved. In R, I know I can create a blank dataframe with 10 lines, then save each read count in a designated row/column. How do I do something similar with a bash loop?

[kubu4]

Without adding any "fancy" stuff like the corresponding filename, you need to change your redirect symbol from > to >>.

The >> will concatenate the results from each run of the loop to the end of the file. Using > will overwrite the file each time, thus, you only end up with the results from the last FastQ file from the loop.

[yaaminiv]

(Posted this in the wrong issue initially...posting here)

Comparing Bismark Processing Reports to the "Total Sequences" from FastQC reports, it looks like they are the same counts...? @kubu4 am I looking at the wrong thing to obtain pre-trimming read counts? Or was there minimal trimming?

[kubu4]

Link to FastQC?

[yaaminiv]

http://owl.fish.washington.edu/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/20180411_fastqc_trim_10bp_Cvirginica_MBD/zr2096_1_s1_R1_val_1_fastqc.html

[kubu4]

That's not the FastQC report for the files in the Bismark report linked above.

[yaaminiv]

Whoops! Here you go:

http://owl.fish.washington.edu/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/20180411_fastqc_trim_10bp_Cvirginica_MBD/zr2096_2_s1_R1_val_1_fastqc.html

http://owl.fish.washington.edu/Athaliana/20180411_trimgalore_10bp_Cvirginica_MBD/20180411_fastqc_trim_10bp_Cvirginica_MBD/zr2096_2_s1_R2_val_2_fastqc.html

[kubu4]

Yes, you are looking at the wrong FastQC report. You're looking at the FastQC report for the trimmed files.

Check nightingales spreadsheet for link to raw FastQ files and/or FastQC for raw files (scroll all the way to right on that spreadsheet to see if there's already a link to a FastQC report for a given FastQ file).

[yaaminiv]

Gotcha. I keep forgetting about the nightingales spreadsheet.

To clarify: Read1URL is the original file, and Read1FastQC is the trimmed file? (can't find that information in the Data Management wiki...is that metadata somewhere else that I don't know of?)