Skip to content
This repository has been archived by the owner on Jan 28, 2023. It is now read-only.

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@CTLife
CTLife authored Jul 5, 2020
1 parent a6a3462 commit a95dca9
Showing 1 changed file with 383 additions and 0 deletions.
383 changes: 383 additions & 0 deletions m6A-SAC-seq_5_SE.pl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,383 @@
#!/usr/bin/env perl
use strict;
use warnings;
use v5.12;
## Perl5 version >= 5.12
###################################################################################################################################################################################################





###################################################################################################################################################################################################
my $input_g = ''; ## such as "4-removedBarcodes"
my $output_g = ''; ## such as "5-finalFASTQ"
{
## Help Infromation
my $HELP = '
------------------------------------------------------------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------------------------------------------------------------
Step 5: Further trim reads by using bbduk.sh of BBMap.
And directly assess the quality of raw reads to identify possible experimental and sequencing errors or biases
by using 4 softwares: FastQC, fastp, FastQ_Screen, and MultiQC.
If this script works well, you do not need to check the the versions of the softwares or packages whcih are used in this script.
And you do not need to exactly match the versions of the softwares or packages.
If some errors or warnings are reported, please check the versions of softwares or packages.
The versions of softwares and packages are used in this script (This is shown for reference only. And do not need to match them.):
Perl, 5.26.1 (perl -v)
FastQC, v0.11.9 (fastqc -v)
fastp, 0.20.0 (fastp -v)
FastQ_Screen, 0.14.0 (fastq_screen -v)
bbduk.sh, 38.76 (bbduk.sh -v)
MultiQC, 1.8 (multiqc --version)
Usage:
m6A-SAC-seq_5_SE.pl [-version] [-help] [-in inputDir] [-out outDir]
For instance:
m6A-SAC-seq_5_SE.pl -in 4-removedBarcodes -out 5-finalFASTQ > m6A-SAC-seq_5_SE.runLog.txt 2>&1
----------------------------------------------------------------------------------------------------------
Optional arguments:
-version Show version number of this program and exit.
-help Show this help message and exit.
Required arguments:
-in inputDir "inputDir" is the name of input path that contains your FASTQ files. (no default)
-out outDir "outDir" is the name of output path that contains your running results (fastq files) of this step. (no default)
-----------------------------------------------------------------------------------------------------------
For more details about this pipeline and other NGS data analysis piplines, please visit https://github.com/CTLife/2ndGS_Pipelines
------------------------------------------------------------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------------------------------------------------------------
';

## Version Infromation
my $version = " The 5th Step, version 1.0.0, 2020-03-01.";

## Keys and Values
if ($#ARGV == -1) { say "\n$HELP\n"; exit 0; } ## when there are no any command argumants.
if ($#ARGV%2 == 0) { @ARGV = (@ARGV, "-help") ; } ## when the number of command argumants is odd.
my %args = @ARGV;

## Initialize Variables
$input_g = '4-removedBarcodes'; ## This is only an initialization value or suggesting value, not default value.
$output_g = '5-finalFASTQ'; ## This is only an initialization value or suggesting value, not default value.

## Available Arguments
my $available = " -version -help -in -out ";
my $boole = 0;
while( my ($key, $value) = each %args ) {
if ( ($key =~ m/^\-/) and ($available !~ m/\s$key\s/) ) {say "\n\tCann't recognize $key"; $boole = 1; }
}
if($boole == 1) {
say "\tThe Command Line Arguments are wrong!";
say "\tPlease see help message by using 'm6A-SAC-seq_5_SE.pl -help' \n";
exit 0;
}

## Get Arguments
if ( exists $args{'-version' } ) { say "\n$version\n"; exit 0; }
if ( exists $args{'-help' } ) { say "\n$HELP\n"; exit 0; }
if ( exists $args{'-in' } ) { $input_g = $args{'-in' }; }else{say "\n -in is required.\n"; say "\n$HELP\n"; exit 0; }
if ( exists $args{'-out' } ) { $output_g = $args{'-out' }; }else{say "\n -out is required.\n"; say "\n$HELP\n"; exit 0; }

## Conditions
$input_g =~ m/^\S+$/ || die "\n\n$HELP\n\n";
$output_g =~ m/^\S+$/ || die "\n\n$HELP\n\n";

## Print Command Arguments to Standard Output
say "\n
################ Arguments ###############################
Input Path: $input_g
Output Path: $output_g
###############################################################
\n";
}
###################################################################################################################################################################################################





###################################################################################################################################################################################################
say "\n\n\n\n\n\n##################################################################################################";
say "Running ......";
sub myMakeDir {
my $path = $_[0];
if ( !( -e $path) ) { system("mkdir -p $path"); }
if ( !( -e $path) ) { mkdir $path || die; }
}
my $output2_g = "$output_g/QC_Results";
&myMakeDir("$output_g");
&myMakeDir("$output2_g");
opendir(my $DH_input_g, $input_g) || die;
my @inputFiles_g = readdir($DH_input_g);
my $numCores_g = 8;
###################################################################################################################################################################################################





###################################################################################################################################################################################################

## These commands must be available:
say "\n\n\n\n\n\n##################################################################################################";
say "Checking all the necessary softwares in this step ......";

sub printVersion {
my $software = $_[0];
system("echo '##############################################################################' >> $output2_g/VersionsOfSoftwares.txt 2>&1");
system("echo '#########$software : ' >> $output2_g/VersionsOfSoftwares.txt 2>&1");
system("$software >> $output2_g/VersionsOfSoftwares.txt 2>&1");
system("echo '\n\n\n\n\n\n' >> $output2_g/VersionsOfSoftwares.txt 2>&1");
}
sub fullPathApp {
my $software = $_[0];
say($software);
system("which $software > yp_my_temp_1.$software.txt");
open(tempFH, "<", "yp_my_temp_1.$software.txt") or die;
my @fullPath1 = <tempFH>;
($#fullPath1 == 0) or die;
system("rm yp_my_temp_1.$software.txt");
$fullPath1[0] =~ s/\n$// or die;
return($fullPath1[0]);
}

&printVersion(" bbduk.sh -v ");
&printVersion(" fastqc --version ");
&printVersion(" fastp --version ");
&printVersion(" fastq_screen --version ");
&printVersion(" multiqc --version ");

my $Trimmomatic_g = &fullPathApp("trimmomatic.jar");
$Trimmomatic_g =~ s/trimmomatic.jar$/Illumina_Small_RNA_5_Adapter.fasta/ or die;
###################################################################################################################################################################################################





###################################################################################################################################################################################################
say "\n\n\n\n\n\n##################################################################################################";
say "Detecting single-end and paired-end FASTQ files ......";
my @singleEnd_g = ();
my @pairedEnd_g = ();
open(seqFiles_FH_g, ">", "$output2_g/singleEnd-pairedEnd-Files.txt") or die;

for ( my $i=0; $i<=$#inputFiles_g; $i++ ) {
next unless $inputFiles_g[$i] !~ m/^[.]/;
next unless $inputFiles_g[$i] !~ m/[~]$/;
next unless $inputFiles_g[$i] !~ m/\.sh$/;
next unless $inputFiles_g[$i] !~ m/^unpaired/;
next unless $inputFiles_g[$i] !~ m/^QC_Results$/;
next unless $inputFiles_g[$i] =~ m/\.fastq/;
my $temp = $inputFiles_g[$i];
say "\t......$temp";
if ($temp !~ m/^(\S+)\.R[12]\.fastq/) { ## sinlge end sequencing files.
$temp =~ m/^(\S+)\.fastq/ or die;
$singleEnd_g[$#singleEnd_g+1] = $temp;
say "\t\t\t\tSingle-end sequencing files: $temp\n";
say seqFiles_FH_g "Single-end sequencing files: $temp\n";
}else{ ## paired end sequencing files.
$temp =~ m/^(\S+)\.R[12]\.fastq/ or die;
$pairedEnd_g[$#pairedEnd_g+1] = $temp;
say "\t\t\t\tPaired-end sequencing files: $temp\n";
say seqFiles_FH_g "Paired-end sequencing files: $temp\n";
}
}

@singleEnd_g = sort(@singleEnd_g);
@pairedEnd_g = sort(@pairedEnd_g);

( ($#pairedEnd_g+1)%2 == 0 ) or die;
say seqFiles_FH_g "\n\n\n\n\n";
say seqFiles_FH_g "All single-end sequencing files:@singleEnd_g\n\n\n";
say seqFiles_FH_g "All paired-end sequencing files:@pairedEnd_g\n\n\n";
say "\t\t\t\tAll single-end sequencing files:@singleEnd_g\n\n";
say "\t\t\t\tAll paired-end sequencing files:@pairedEnd_g\n\n";

my $numSingle_g = $#singleEnd_g + 1;
my $numPaired_g = $#pairedEnd_g + 1;

say seqFiles_FH_g "\nThere are $numSingle_g single-end sequencing files.\n";
say seqFiles_FH_g "\nThere are $numPaired_g paired-end sequencing files.\n";
say "\t\t\t\tThere are $numSingle_g single-end sequencing files.\n";
say "\t\t\t\tThere are $numPaired_g paired-end sequencing files.\n";


for ( my $i=0; $i<$#pairedEnd_g; $i=$i+2 ) {
my $temp = $pairedEnd_g[$i];
$temp =~ s/\.R1\.fastq/.R2.fastq/ or die "\n##Error-1: $temp ##\n\n";
($pairedEnd_g[$i+1] eq $temp) or die "\n##Error-2: $temp ## $pairedEnd_g[$i+1] ##\n\n";
}

print("\n\n\n\n\n#########################################\n");
###################################################################################################################################################################################################





###################################################################################################################################################################################################
{
say "\n\n\n\n\n\n##################################################################################################";
say "Trim reads by using bbduk.sh of BBMap ......";
for (my $i=0; $i<=$#pairedEnd_g; $i=$i+2) {
$pairedEnd_g[$i] =~ m/\.R1.fastq/ or die;
my $temp = $pairedEnd_g[$i];
my $end1 = $temp;
my $end2 = $temp;
$end2 =~ s/\.R1\.fastq/.R2.fastq/ or die "\n##Error-3: $temp ##\n\n";
say "\t......$end1";
say "\t......$end2";
($end2 eq $pairedEnd_g[$i+1]) or die;
open(tempFH, ">>", "$output2_g/paired-end-files.txt") or die;
say tempFH "$end1, $end2\n";
$temp =~ s/\.fastq\.\S+$// ;
$temp =~ s/\.fastq$// ;
system("bbduk.sh in=$input_g/$end1 in2=$input_g/$end2 out=$output_g/$end1 out2=$output_g/$end2 ref=$Trimmomatic_g stats=$output2_g/$temp.stats statscolumns=5 minlength=10 ktrim=l k=12 hdist=0 >> $output2_g/$temp.runLog 2>&1"); ##mink=11
}
for (my $i=0; $i<=$#singleEnd_g; $i++) {
say "\t......$singleEnd_g[$i]" ;
my $temp = $singleEnd_g[$i];
say "\t......$temp";
$temp =~ s/\.fastq\.\S+$// ;
$temp =~ s/\.fastq$// ;
system("bbduk.sh in=$input_g/$singleEnd_g[$i] out=$output_g/$singleEnd_g[$i] ref=$Trimmomatic_g stats=$output2_g/$temp.stats statscolumns=5 minlength=12 ktrim=l k=12 hdist=0 >> $output2_g/$temp.runLog.txt 2>&1"); ##mink=11
}
}
###################################################################################################################################################################################################





###################################################################################################################################################################################################
sub myQC_FASTQ_1 {
my $dir1 = $_[0]; ## All the fastq files must be in this folder.
my $QCresults = "$dir1/QC_Results";
my $FastQC = "$QCresults/1_FastQC";
my $fastp = "$QCresults/2_fastp";
my $MultiQC1 = "$QCresults/MultiQC/1_FastQC";
my $MultiQC2 = "$QCresults/MultiQC/2_fastp";
&myMakeDir($QCresults);
&myMakeDir($FastQC);
&myMakeDir($fastp);
&myMakeDir($MultiQC1);
&myMakeDir($MultiQC2);
opendir(my $FH_Files, $dir1) || die;
my @Files = readdir($FH_Files);
say "\n\n\n\n\n\n##################################################################################################";
say "Detecting the quality of all FASTQ files by using FastQC, fastp and MultiQC ......";
for ( my $i=0; $i<=$#Files; $i++ ) {
next unless $Files[$i] =~ m/\.fastq/;
next unless $Files[$i] !~ m/^[.]/;
next unless $Files[$i] !~ m/[~]$/;
next unless $Files[$i] !~ m/\.sh$/;
next unless $Files[$i] !~ m/^unpaired/;
next unless $Files[$i] !~ m/^QC_Results$/;
my $temp = $Files[$i];
say "\t......$temp";
$temp =~ s/\.fastq\.\S+$// ;
$temp =~ s/\.fastq$// ;
system( "fastqc --outdir $FastQC --threads $numCores_g --format fastq --kmers 6 $dir1/$Files[$i] >> $FastQC/$temp.runLog.txt 2>&1" );
system( "fastp --in1 $dir1/$Files[$i] -p --report_title $temp --thread $numCores_g --json $fastp/$temp.fastp.json --html $fastp/$temp.fastp.html --report_title $temp --reads_to_process 100000 >> $fastp/$temp.runLog.txt 2>&1" );
}
system( "multiqc --title FastQC --filename FastQC --module fastqc --verbose --export --outdir $MultiQC1 $FastQC >> $MultiQC1/MultiQC.FastQC.runLog.txt 2>&1" );
system( "multiqc --title fastp --filename fastp --module fastp --verbose --export --outdir $MultiQC2 $fastp >> $MultiQC2/MultiQC.fastp.runLog.txt 2>&1" );
}
###################################################################################################################################################################################################





###################################################################################################################################################################################################
sub myQC_FASTQ_2 {
my $dir1 = $_[0]; ## All the fastq files must be in this folder.
my $QCresults = "$dir1/QC_Results";
my $FastQ_Screen = "$QCresults/3_FastQ-Screen";
my $MultiQC = "$QCresults/MultiQC/3_FastQ-Screen";
&myMakeDir($QCresults);
&myMakeDir($FastQ_Screen);
&myMakeDir($MultiQC);
opendir(my $FH_Files, $dir1) || die;
my @Files = readdir($FH_Files);
say "\n\n\n\n\n\n##################################################################################################";
say "Detecting the quality of all FASTQ files by using FastQ_Screen and MultiQC ......";
for ( my $i=0; $i<=$#Files; $i++ ) {
next unless $Files[$i] =~ m/\.fastq/;
next unless $Files[$i] !~ m/^[.]/;
next unless $Files[$i] !~ m/[~]$/;
next unless $Files[$i] !~ m/\.sh$/;
next unless $Files[$i] !~ m/^unpaired/;
next unless $Files[$i] !~ m/^QC_Results$/;
my $temp = $Files[$i];
say "\t......$temp";
$temp =~ s/\.fastq\.\S+$// ;
$temp =~ s/\.fastq$// ;
system( "fastq_screen --aligner bowtie2 --outdir $FastQ_Screen/$temp --threads $numCores_g --subset 100000 $dir1/$Files[$i] >> $FastQ_Screen/$temp.runLog.txt 2>&1" );
}
system( "multiqc --title FastQ_Screen --filename FastQ_Screen --module fastq_screen --verbose --export --outdir $MultiQC $FastQ_Screen >> $MultiQC/MultiQC.FastQ-Screen.runLog.txt 2>&1" );
}
###################################################################################################################################################################################################





###################################################################################################################################################################################################
&myQC_FASTQ_1($output_g);
&myQC_FASTQ_2($output_g);
###################################################################################################################################################################################################





###################################################################################################################################################################################################
{
say "\n\n\n\n\n\n##################################################################################################";
say "Detecting the quality of paired-end FASTQ files by using fastp and FastQ_Screen with PE mode ......";
my $fastp = "$output2_g/4_fastp_PairedEnd";
my $MultiQC1 = "$output2_g/MultiQC/4_fastp_PairedEnd";
&myMakeDir($fastp);
&myMakeDir($MultiQC1);
for ( my $j=0; $j<=$#pairedEnd_g; $j=$j+2 ) {
my $temp1 = $pairedEnd_g[$j];
my $temp2 = $pairedEnd_g[$j+1];
my $temp = $temp2;
say "\t......$temp1";
say "\t......$temp2";
$temp1 =~ m/\.R1\.fastq/ or die;
$temp2 =~ m/\.R2\.fastq/ or die;
$temp =~ s/\.R2\.fastq\.\S+$// ;
$temp =~ s/\.R2\.fastq$// ;
system( "fastp --in1 $output_g/$temp1 --in2 $output_g/$temp2 -p --report_title $temp --thread $numCores_g --json $fastp/$temp.fastp.json --html $fastp/$temp.fastp.html >> $fastp/$temp.runLog 2>&1" );
}
system( "multiqc --title fastp --filename fastp --module fastp --verbose --export --outdir $MultiQC1 $fastp >> $MultiQC1/MultiQC.fastp.runLog.txt 2>&1" );
}
###################################################################################################################################################################################################





###################################################################################################################################################################################################
say "\n\n\n\n\n\n##################################################################################################";
say "\tJob Done! Cheers! \n\n\n\n\n";





## END




0 comments on commit a95dca9

Please sign in to comment.