Pipeline successfully ran:
Unable to render embedded object: File (Screen Shot 2023-02-09 at 9.40.46 AM.png) not found.

Directory: /nobackup/tomato_genome/30-804059537-KP

Comment: There are no errors in the report but the number of sequences mapped is pretty low. Might need to look into that! Double check sample sheet I made maybe or the wrong reference genome was used to map data.

• Sequence Duplication levels might be an issue
  Unable to render embedded object: File (Screen Shot 2023-02-09 at 10.28.15 AM.png) not found.
• Per sequence GC Content also poor
• Alignment scores are poor because unmapped reads are too short.
  
  Link to interpret report: https://nf-co.re/eager/2.2.2/output#multiqc-report

Next steps:

• remove sorted names
• make coverage graphs

The "star alignment scores" section looks informative. According to the plot, a high percentage of sequences in some samples were reported as "Unmapped: too short." I think this can happen when the library contained a lot of very short inserts. I think we should go ahead and proceed with the rest of the pipeline, but keep an eye on how those samples with the largest percentage of "shorties" perform in subsequent analyses.

attn: [~molly]

Other people were having the same issue and said:

"I encountered this issue when in the two paired-end input FASTQ files mates are out-of-order, i.e. mates are not found at the same line of the two files. This leads to a lot of not properly mapped read pairs that STAR throws into the "too short" bucket." https://github.com/alexdobin/STAR/issues/169

Not sure if this is the same but I will look more into it and check the fastq file directory. The naming of the files might be confusing nextflow as well because R1 and R2 are used twice in some file names.

[~aloraine]

Renaming code: https://bitbucket.org/hotpollen/splicing-analysis/src/main/src/renameBams.sh

Update:

• I have changed the names of the fastq files so instead of having R1, R2, R3 it is now Rep1, Rep2, Rep3. This is due to nextflow possibly confusing paired end files and not being matched together correctly.
• Here is the new csv samples file: [^KP_samples.csv]
• Rerunning nextflow.

Update:

New MutliQC report:
[^KP_multiqc_report.html]

Comment: Unfortunately, even after changing the fastq file names the mutliqc report is the exact same as the last one with Unmapped: too short # of reads for the alignment scores.

Here is the txt file to see the alignment scores instead of using the actual graph from the report:
multiqc_star.txt

Directory: /nobackup/tomato_genome/30-804059537-KP/results/multiqc/star_salmon/multiqc_data

Next steps:

• Proceed with using the output from most recent run of nextflow nf-core/rnaseq pipeline which used the renamed samples (e.g., Rep1, Rep2, Rep3)
• [~aloraine] to review samples file for possible problems

Update:

I moved forward with results and renamed sorted bam files and made coverage graphs.

Directory: /nobackup/tomato_genome/30-804059537-KP/results/star_salmon

Update: After speaking with Nowlan, we agree that running fastqc on the files would be beneficial to check the quality of the data before running it with nextflow.

Script:

#!/bin/bash


#SBATCH --job-name=fastqc               #job name after submission
#SBATCH -p Orion                        #partition being used
#SBATCH -N 1                            #number of nodes to use
#SBATCH --ntasks-per-node=8             #max number of tasks per node
#SBATCH --mem=60gb                      #memory required per node
#SBATCH -t 0-50:00                      #time (D-HH:MM)
#SBATCH -o fastqc.%j.out                #standard output file
#SBATCH -e fastqc.%j.err                #standard error file
#SBATCH --mail-type=END,FAIL            #Notifications for job complete/failure
#SBATCH --mail-user=mdavi258@uncc.edu   #Send to user email


module load fastqc

for i in /nobackup/tomato_genome/30-804059537-KP/*.fastq.gz
do
  	fastqc -o /nobackup/tomato_genome/30-804059537-KP/fastQC_dataQuality $i
done

Directory: /nobackup/tomato_genome/30-804059537-KP/fastQC_dataQuality

For quick reference here are two of the fastqc reports to check data quality:
[^Heinz-Ovary-Rep1-0hr-25C-unpol_R1_001_fastqc.html]
[^Heinz-Ovary-Rep2-0hr-25C-unpol_R1_001_fastqc.html]

Let me know what you think Nowlan Freese

Notes about High Duplication:

• High duplication is either going to be the result of technical duplication (too many PCR cycles), or over-sequencing (very high fold coverage).
• How many PCR cycles were done during the protocol? Also, how many reads were their total?
• Dimer Contamination?
• Should we just let it through for downstream analysis?

Source: https://wiki.bits.vib.be/index.php/Quality_control_of_NGS_data

Validating mate pair files

Script:

#!/bin/bash


#SBATCH --job-name=validate_script      #job name after submission
#SBATCH -p Orion                        #partition being used
#SBATCH -N 1                            #number of nodes to use
#SBATCH --ntasks-per-node=8             #max number of tasks per node
#SBATCH --mem=60gb                      #memory required per node
#SBATCH -t 0-50:00                      #time (D-HH:MM)
#SBATCH -o validate.%j.out              #standard output file
#SBATCH -e validate.%j.err              #standard error file
#SBATCH --mail-type=END,FAIL            #Notifications for job complete/failure
#SBATCH --mail-user=mdavi258@uncc.edu   #Send to user email
#SBATCH --array=1-63


#setting up where to grab files from
file=$(sed -n -e "${SLURM_ARRAY_TASK_ID}p"  /nobackup/tomato_genome/30-804059537-KP/kp_runlist.txt)

#The command to validate each pair:
cd /nobackup/tomato_genome/30-804059537-KP

perl /projects/tomato_genome/scripts/validateHiseqPairs.pl ${file}_R1_001.fastq.gz ${file}_R2_001.fastq.gz


echo "Done"

Error: mate-pair files are not ordered

Ann's note: We need to understand what exactly this script actually is doing and assessing.

[~RobertReid] : MD5 checking shows that the fastq files were not corrupted by the transfer.
Nowlan Freese : Suggests comparing the first few lines per file (via "head" function) to see if pairs are present

Problem is: The mate pair records in the "1" and "2" fastq files per sample are not in the same order in the two files. According to the above script (validating mate pairs) the records are out of order.

MD5 Check:

All the MD5 check out.

Details on this are saved to Google Drive in Kelsie's experiment folder.
https://drive.google.com/drive/folders/1TxUDhJHr9mrXOVysrcceS9YGyTXFHRmo?usp=share_link

Please add validateMatePairs.pl to "src" directory in the repository:

• https://bitbucket.org/hotpollen/pistil-rna-seq/src/main/.

NF suggestion: Run fastqc on output of trim galore to observe new and possibly aberrant size distribution of read sequence

validateMatePairs.pl code:

#!/usr/bin/perl

use strict;
use warnings;

open FH1,$ARGV[0] or die "\n can not open file $ARGV[0]\n";  ## first pair
open FH2,$ARGV[1] or die "\n can not open file $ARGV[1]\n";  ## second pair
my($str1,$str2,$tempStr);
my($n1,$n2);
$n1 = 0;
$n2 = 0;
print " Validating $ARGV[0] and $ARGV[1] \n";
my(@a1,@a2);

	while($str1 = <FH1>){
        $tempStr = <FH1>;
        $tempStr = <FH1>;
        $tempStr = <FH1>;
        ++$n1;

        $str2 = <FH2>;
        $tempStr = <FH2>;
        $tempStr = <FH2>;
        $tempStr = <FH2>;
        ++$n2;

        $str1 =~ s/\n//;
        $str1 =~ s/\r//;
        $str2 =~ s/\n//;
        $str2 =~ s/\r//;


        @a1 = split(/\s+/,$str1);
        @a2 = split(/\s+/,$str2);

        $str1 = $a1[0];
        $str2 = $a2[0];

        $str1 =~ s/(\/\d)$//;
        $str2 =~ s/(\/\d)$//;

                if($str1 ne $str2){
                die "Read pairs not found for $str1 and $str2, mate-pair files are not ordered\n";
                }
        }  ## while(<FH1>) ends
close FH1;
close FH2;

print "Total validated mates: $n1 and $n2: Read-pairs are properly ordered\n";

Next step: Run the code with unzipped fastq files.

gzip -d *.gz

Comment: script finished running and output files say Read-pairs are properly ordered. Decompressing the files helped fix the validate script error.

I ran trimmomatic on the raw data to see what would happen.

Resulting files are located here:
/nobackup/tomato_genome/30-804059537-KP/trimmoTest

Script to run it is here:
/projects/tomato_genome/scripts/rob/trimmomatic-temp.slurm

The script will validate pairs afterwards. (using option -validate pairs)
It appears that 99.95% of all reads are great.
About 2-8 reads are bad per pairing.
Those reads are saved as unpaired.fastq files.

We could now run nextflow with these read files instead.

[~RobertReid] suggests counting number of aligned versus unaligned reads, for each BAM file we have.

To use samtools to view the aligned and unaligned.

READS MAPPED:
module load samtools
samtools view -c -F 4 nagcarlang-sorted.bam

For Unmapped:
samtools view -c -f 4 nagcarlang-sorted.bam

Let's calculate coverage:
samtools depth nagcarlang-sorted.bam | awk '

END

'

Put these lines into a slurm script. Should run very quickly.

Created pull request to add src folder and perl validation file to bitbucket:

https://bitbucket.org/hotpollen/pistil-rna-seq/pull-requests/1

[~aloraine]

Created script to use samtools to view the aligned and unaligned:

#!/bin/bash


#SBATCH --job-name=samtools_view        #job name after submission
#SBATCH -p Orion                        #partition being used
#SBATCH -N 1                            #number of nodes to use
#SBATCH --ntasks-per-node=8             #max number of tasks per node
#SBATCH --mem=900gb                     #memory required per node
#SBATCH -t 14-00:00                     #time (D-HH:MM)
#SBATCH -o samtools_view.%j.out         #standard output file
#SBATCH --mail-type=END,FAIL            #Notifications for job complete/failure
#SBATCH --mail-user=mdavi258@uncc.edu   #Send to user email
#SBATCH --array=1-63


file=$(sed -n -e "${SLURM_ARRAY_TASK_ID}p"  /nobackup/tomato_genome/30-804059537-KP/kp_runlist.txt)

module load samtools
echo "Mapped:" ${file}
samtools view -c -F 4 ${file}.bam
echo  "Unmapped:" ${file}
samtools view -c -f 4 ${file}.bam

echo "Calculate Coverage" 
samtools depth ${file}.bam | awk '{sum+=$3} END { print "Average = ",sum/NR}'

echo "done"
echo "---------------------------------------------------------"

Directory: /nobackup/tomato_genome/30-804059537-KP/results/star_salmon

Combined output files into one:

cat *.out > ./mergedsamtoolsOut.txt

Output File:

mergedsamtoolsOut.txt

Ann's comments:

Based on output above:

• The bam files do not contain any unmapped reads, only mapped reads
• The samtools "depth" command computes the number of alignments per base pair position - see http://www.htslib.org/doc/samtools-depth.html
• For transcriptome data, the "depth" command does not make a lot of sense because the depth of read alignments at any given position depends on whether or not that position is inside an exon, and also on the level of expression of that exon
• I don't know what "NR" means and where this is coming from in the "sum/NR" statement at the end of the script

Conclusion: This output of this script script does not explain the QC result.

We do not know why some of the samples did not perform well. Let's proceed with the pipeline and visualize the data in a genome browser as this visualization step may reveal more information about the problematic samples.

I am including the Google Links that Kelsie provided here for documentation purposes.

Tomato Pistil Tissue Collection Protocol
https://docs.google.com/document/d/1g8GJBEzxUC-QjfMXk0Eq5mT8e31bGirjXYM3-Sv4u7Q/edit?usp=sharing

Experimental Design for Solavar
https://docs.google.com/document/d/1BXVq-0oop3Ch3Qzbr2nkhQczG1cZH2yBqQRMmBXKGyo/edit?usp=sharing

Sequenced Samples
https://docs.google.com/spreadsheets/d/1WwPzifPzbACmgS3uR_V92cYIGN3qS1yWPGHBu7DY_-I/edit?usp=sharing