###RNA seq data processing

subread-buildindex -o /PATH/Indexes/At_Genome_subread /PATH/Genomes/FASTA/TAIR10.fasta

for PSET in SG335 SG337 SG339B SG339C SG339E SG340A SG340B SG342 SG346 SG350 SG352 SG353 SG354 SG355 SG356 SG357 SG358 SG359 SG362 SG366 SG367 SG369 SG373
do mkdir /PATH/RNAseq_alignment/"$PSET"
done


for PSET in SG335 SG337 SG339B SG339C SG339E SG340A SG340B SG342 SG346 SG350 SG352 SG353 SG354 SG355 SG356 SG357 SG358 SG359 SG362 SG366 SG367 SG369 SG373 
do echo "$PSET"
subread-align -i /PATH/Genomes/Indexes/At_Genome_subread -r /PATH/RNAseqExtracted/"$PSET"_extracted/"$PSET"merged.fastq --BAMoutput -o /PATH/RNAseq_alignment/"$PSET"/"$PSET"_aligned.bam
done

##sorting and indexing
for PSET in SG335 SG337 SG339B SG339C SG339E SG340A SG340B SG342 SG346 SG350 SG352 SG353 SG354 SG355 SG356 SG357 SG358 SG359 SG362 SG366 SG367 SG369 SG373 
do echo "$PSET"
samtools sort -m 12000000000 /PATH/RNAseq_alignment/"$PSET"/"$PSET"_aligned.bam /PATH/RNAseq_alignment/"$PSET"/"$PSET"_aligned_sorted
done
for PSET in SG335 SG337 SG339B SG339C SG339E SG340A SG340B SG342 SG346 SG350 SG352 SG353 SG354 SG355 SG356 SG357 SG358 SG359 SG362 SG366 SG367 SG369 SG373 
do echo "$PSET"
samtools index /PATH/RNAseq_alignment/"$PSET"/"$PSET"_aligned_sorted.bam
done


##use Rcount
RcountCLVPATH=/home/stefan/ownCloud/Bioinformatik/Rcount/linux_64bit/Rcount_CLV
for PSET in SG339B SG339C SG339E SG340A SG340B SG342 SG346 SG350 SG352 SG353 SG354 SG355 SG356 SG357 SG358 SG359 SG362 SG366 SG367 SG369 SG373
do echo "$PSET"
"$RcountCLVPATH"/Rcount-multireads_CLV -c /PATH/RNAseq_alignment/"$PSET"/"$PSET"_aligned_sorted.bam,/PATH/RNAseq_alignment/"$PSET"/"$PSET"_aligned_multiout.bam,y,100
done


#### get xml project files

RcountCLVPATH="/PATH/Rcount/linux_64bit/Rcount_CLV"
RCOUNTXML="/PATH/Rcount_AT_pROK_detailed.xml"
multioutPATH="/PATH/alignedTo_genome_with_pROK"
OUTpath="/PATH/MappedToAtGenome_pROK/RcountDistribute_OUT"
for PSET in SG335 SG337 SG339B SG339C SG339E SG340A SG340B SG342 SG346 SG350 SG352 SG353 SG354 SG355 SG356 SG357 SG358 SG359 SG362 SG366 SG367 SG369 SG373 
do echo $PSET
python $RcountCLVPATH/prepareXMLfileForRcount.py $RCOUNTXML $multioutPATH/"$PSET"/"$PSET"_aligned_multiout.bam $OUTpath/"$PSET"_multiout_mapped.bam $OUTpath/"$PSET"_Rcount_counts.txt "$PSET" true true 5 250 1 > $OUTpath/"$PSET"_Rcount_project.xml
done

for PSET in SG335 SG337 SG339B SG339C SG339E SG340A SG340B SG342 SG346 SG350 SG352 SG353 SG354 SG355 SG356 SG357 SG358 SG359 SG362 SG366 SG367 SG369 SG373
do echo $PSET
$RcountCLVPATH/Rcount-distribute-antisense -c $OUTpath/"$PSET"_Rcount_project.xml
done


##########################################################################################################################################################################
############################################################# smRNAseq ###################################################################################################
##########################################################################################################################################################################
FASTQPATH=/PATH/smRNA
trimmedFASTQ=/PATH/trimmedFastq
BOWTIEPATH_pROK2=/PATH/Indexes/BOWTIEpROK2
RESPATH=/PATH/aligned_sRNA
RcountCLVPATH=/PATH/Rcount/linux_64bit/Rcount_CLV
ANALYSISPATH=/PATH/smRNAseq/analysis
PYTPATH=/PATH/

##cut adapters
for seq in 261 335 337 339A 339B 339C 340 342 346 350 352 353 354 355 356 357 358 359 361 362 368 369 371
do echo "$seq"
cutadapt -b TGGAATTCTCGGGTGCCAAGGAACTCCAGTCAC -m 17 -q 20 "$FASTQPATH"/20180615.A-"$seq"_R1.fastq -o "$trimmedFASTQ"/SG"$seq"_smRNA_trimmed.fastq
done


grep -c TGGAATTCTCGGGTGCCAAGGAACTCCAGTCAC "$FASTQPATH"/20180615.A-261_R1.fastq

grep  -c TCGCTTGGTGCAGGTCGGGAA "$FASTQPATH"/20180503.A-261_R1.fastq
wc -l "$FASTQPATH"/20180503.A-261_R1.fastq

grep  -m 100 TGGAATTCTCGGGTGCCAAGGAACTCCAGTCAC "$FASTQPATH"/20180503.A-261_R1.fastq


#here we map against rRNA and tRNA
#fasta files for all rRNA and tRNA loci where made in R using system(python...)
#we add the two weird regions to the fasta...
python PYTPATH/FASTAseqExtract.py Chr2 1 10000 /PATH/TAIR10_allChromosomes_noNewLine.fasta /PATH/rRNA_artefact_Chr2.fasta
python PYTPATH/FASTAseqExtract.py Chr3 14194000 14204000 /PATH/TAIR10_allChromosomes_noNewLine.fasta /PATH/rRNA_artefact_Chr3.fasta
cat /media/DATA3/Genomes/TAIR10/rRNA_tRNA_TAIR10.fasta /PATH/rRNA_artefact_Chr2.fasta /PATH/rRNA_artefact_Chr3.fasta /PATH/AT_organelle.fas > /PATH/regionsToFilter.fasta

bowtie-build /media/DATA3/Genomes/TAIR10/regionsToFilter.fasta /media/DATA3/Genomes/Indexes/ToFilter/regionsToFilter
#now I filter out reads that map to rRNA tRNA and the weird regions after rRNAs on chrom3 and chrom2 and the chloroplasts and mitochondria
for seq in 261 335 337 339A 339B 339C 340 342 346 350 352 353 354 355 356 357 358 359 361 362 368 369 371
do echo SG"$seq"
bowtie -v 2 --best -p 20 -m 1000 -t /PATH/regionsToFilter -q "$trimmedFASTQ"/SG"$seq"_smRNA_trimmed.fastq --un "$trimmedFASTQ"/SG"$seq"_smRNA_trimmed_clean.fastq
done

##filter for length
for seq in 261 335 337 339A 339B 339C 340 342 346 350 352 353 354 355 356 357 358 359 361 362 368 369 371
do awk 'BEGIN {OFS = "\n"} {header = $0 ; getline seq ; getline qheader ; getline qseq ; if (length(seq) >= 17 && length(seq) <= 30) {print header, seq, qheader, qseq}}' < "$trimmedFASTQ"/SG"$seq"_smRNA_trimmed_clean.fastq > "$trimmedFASTQ"/SG"$seq"_smRNA_trimmed_cleaner.fastq
done

###map using bowtie
for seq in 261 335 337 339A 339B 339C 340 342 346 350 352 353 354 355 356 357 358 359 361 362 368 369 371
do echo "$seq"
bowtie -v 2 --best -p 20 -m 10000 --sam -t "$BOWTIEPATH_pROK2/TAIR10_pROK2" -q "$trimmedFASTQ"/SG"$seq"_smRNA_trimmed_cleaner_pool.fastq "$RESPATH"/SG"$seq"_aligned.sam >> "$ANALYSISPATH"/SG"$seq"_bowtie_stats.txt 2>&1
python "$PYTPATH/stefanMultiSam.py" "$RESPATH"/SG"$seq"_aligned.sam "$RESPATH"/SG"$seq"_aligned.bam
rm "$RESPATH"/SG"$seq"_aligned.sam
samtools sort -m 12000000000 --threads 4 "$RESPATH"/SG"$seq"_aligned.bam -o "$RESPATH"/SG"$seq"_aligned_sorted.bam
rm "$RESPATH"/SG"$seq"_aligned.bam
samtools index "$RESPATH"/SG"$seq"_aligned_sorted.bam
"$RcountCLVPATH/Rcount-multireads_CLV" -c "$RESPATH"/SG"$seq"_aligned_sorted.bam,"$RESPATH"/SG"$seq"_aligned_multiout.bam,y,25
samtools index "$RESPATH"/SG"$seq"_aligned_multiout.bam
samtools idxstats "$RESPATH"/SG"$seq"_aligned_multiout.bam | awk '{print $1"\t"$3}' > "$ANALYSISPATH"/SG"$seq"_stats_multi.txt
samtools idxstats "$RESPATH"/SG"$seq"_aligned_sorted.bam | awk '{print $1"\t"$3}' > "$ANALYSISPATH"/SG"$seq"_stats.txt
done

for seq in 261 335 337 339A 339B 339C 340 342 346 350 352 353 354 355 356 357 358 359 361 362 368 369 371
do samtools view -h "$RESPATH"/SG"$seq"_aligned_multiout.bam | awk 'length($10) == 21 || length($10) == 24 || $1 ~ /^@/' | samtools view -bS - > "$RESPATH"/SG"$seq"_aligned_multiout_21_24.bam
done

for seq in 261 335 337 339A 339B 339C 340 342 346 350 352 353 354 355 356 357 358 359 361 362 368 369 371
do samtools index "$RESPATH"/SG"$seq"_aligned_multiout_21_24.bam
echo "$seq"
samtools idxstats "$RESPATH"/SG"$seq"_aligned_multiout_21_24.bam  | awk '{print $3}'
done

#count initial reads...
for seq in 261 335 337 339A 339B 339C 340 342 346 350 352 353 354 355 356 357 358 359 361 362 368 369 371
do echo "$seq"
wc -l "$trimmedFASTQ"/SG"$seq"_smRNA_trimmed.fastq
done
##now bin the reads in HiCdat...



