#該流程用于對(duì)SRR數(shù)據(jù)直接分析
#若使用fastq文件分析，則去除前面的對(duì)于SRR轉(zhuǎn)換的過(guò)程

#此位置需要更改為自己的基因組信息
REF="/mnt/f/test/ref/B73_v4"
GTF="/mnt/f/test/ref/Zea_mays.B73_RefGen_v4.47.gtf"
threads=22


mkdir 00Rawdata
#下載原始數(shù)據(jù)  #axel為可以多線程下載，比wget下載速度快一些

axel -n 22 -o ./00Rawdata ftp.sra.ebi.ac.uk/vol1/fastq/SRR603/007/SRR6039707/SRR6039707_1.fastq.gz
axel -n 22 -o ./00Rawdata ftp.sra.ebi.ac.uk/vol1/fastq/SRR603/007/SRR6039707/SRR6039707_2.fastq.gz

#利用md5sum 驗(yàn)證數(shù)據(jù)完整性
#md5sum -c ./00Rawdata/file.md5new.txt > ./00Rawdata/md5check.txt
mkdir 01Cleandata 02Alignment 03Result
mkdir ballgown

#替換成自己的文件
for i in  SRR6039725  SRR6039724
do 
    #轉(zhuǎn)化SRR為fastq文件
    #可以使用多線程
    fasterq-dump --split-3 --threads ${threads} ${i} -O ./00Rawdata/
    #使用單線程
    #fastq-dump --split-3 ${i} -O ./00Rawdata/

    #壓縮文件
    #單線程
    #gzip ./00Rawdata/${i}_1.fastq
    #gzip ./00Rawdata/${i}_2.fastq
    
    #多線程,默認(rèn)使用全部線程
    pigz ./00Rawdata/${i}_1.fastq
    pigz ./00Rawdata/${i}_2.fastq


    # 去接頭
    # trimmomatic 可執(zhí)行文件的位置
    # 單端或雙端測(cè)序，單端為SE， 雙端為PE
    # phred33 或者-phred64 指定堿基質(zhì)量編碼，如果未指定，則自動(dòng)識(shí)別。Illumina 使用-phred33
    # -threads 指定使用的線程數(shù)
    # 指定處理后的數(shù)據(jù)所放置的位置，單端測(cè)序如前所示，雙端測(cè)序如下：  ./00Rawdata/${i}_1.fastq.gz ./00Rawdata/${i}_2.fastq.gz ./01Cleandata/${i}_1_clean.fastq.gz ./01Cleandata/${i}_1_unpaired.fastq.gz ./01Cleandata/${i}_2_clean.fastq.gz ./01Cleandata/${i}_2_unpaired.fastq.gz 
    #LEADING:20 \ 從起始端開始去除低質(zhì)量的堿基，只要一個(gè)堿基低于閾值，就會(huì)切除該堿基，并調(diào)查下一個(gè)堿基
    # TRAILING:20從末端移除低質(zhì)量的堿基，只要一個(gè)堿基質(zhì)量低于與之，就會(huì)切除該堿基，并調(diào)查下一個(gè)堿基
    #SLIDINGWINDOW:4:15  當(dāng)窗口內(nèi)的平均質(zhì)量低于閾值時(shí)，執(zhí)行活動(dòng)窗口提出，通過(guò)考慮多個(gè)堿基，單個(gè)質(zhì)量差的堿基不會(huì)導(dǎo)致刪除了在reads后期高質(zhì)量的數(shù)據(jù)     
    #MINLEN:20 設(shè)置保留reads的最小長(zhǎng)度
    #java -jar /mnt/d/ubuntuSoftware/Trimmomatic-0.39/trimmomatic-0.39.jar PE -phred33 -threads $threads ./00Rawdata/${i}_1.fastq.gz ./00Rawdata/${i}_2.fastq.gz ./01Cleandata/${i}_1_clean.fastq.gz ./01Cleandata/${i}_1_unpaired.fastq.gz ./01Cleandata/${i}_2_clean.fastq.gz ./01Cleandata/${i}_2_unpaired.fastq.gz LEADING:20 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:20

    #fastp使用方法
    fastp --thread ${threads} -i ./00Rawdata/${i}_1.fastq.gz -I ./00Rawdata/${i}_2.fastq.gz  -o ./01Cleandata/${i}_1_clean.fastq.gz -O ./01Cleandata/${i}_2_clean.fastq.gz -j ./01Cleandata/${i}.json -h ./01Cleandata/${i}.html 
    
    #比對(duì)到基因組
    #--dta這個(gè)選項(xiàng)告訴 HISAT2 生成用于基因表達(dá)定量的數(shù)據(jù)（DTA，Documented Transcripts Alignment）,還有專門用于下游cufflink分析的--dta-cufflinks
    #-x 選項(xiàng)后面是參考基因組的索引文件的路徑。在這里，參考基因組的索引文件位于 ../ref/B73_v4 目錄中。
    #-p 選項(xiàng)指定了要使用的線程數(shù)，這里設(shè)置為 22，表示使用 22 個(gè)線程來(lái)執(zhí)行比對(duì)任務(wù)，以提高處理速度。
    #-1,-2 選項(xiàng)用于指定待比對(duì)的雙端測(cè)序數(shù)據(jù)文件的路徑。
    # -S 選項(xiàng)指定了輸出比對(duì)結(jié)果的 SAM 格式文件的路徑和名稱。SAM 文件是一種常見(jiàn)的比對(duì)結(jié)果文件格式，它包含了每個(gè)測(cè)序讀取的比對(duì)信息。
    hisat2 --dta -x ${REF} -p ${threads} -1 ./01Cleandata/${i}_1_clean.fastq.gz -2 ./01Cleandata/${i}_2_clean.fastq.gz -S ./02Alignment/${i}.sam  
    
    #將sam文件轉(zhuǎn)化為bam文件      
    #view 是 samtools 的一個(gè)子命令，用于查看、轉(zhuǎn)換和篩選 SAM/BAM 文件。
    #-@ 指定線程數(shù)
    #-S 選項(xiàng)指定了要處理的輸入文件，即 SAM 格式的比對(duì)結(jié)果文件的路徑。
    # -b 選項(xiàng)告訴 samtools 將輸入文件轉(zhuǎn)換為 BAM 格式。BAM 是一種二進(jìn)制格式，通常比 SAM 文件更緊湊，占用更少的存儲(chǔ)空間，并且更容易處理。
    #> 這部分代碼將轉(zhuǎn)換后的 BAM 文件的輸出路徑和文件名指定為 ./02Alignment/${i}.bam。
    samtools view -@ ${threads} -S ./02Alignment/${i}.sam -b > ./02Alignment/${i}.bam  
    
    #為了節(jié)約空間，可以將不用的sam文件刪除
    sudo rm -r ./02Alignment/${i}.sam
    
    #sort 是 samtools 的一個(gè)子命令，用于對(duì) BAM 文件進(jìn)行排序。
        #-@指定線程數(shù)
        # 這是要排序的輸入文件，即未排序的 BAM 文件的路徑。${i} 變量似乎用于迭代多個(gè)樣本文件，因此 ${i} 會(huì)在循環(huán)中被實(shí)際的文件名替換。這里假設(shè)未排序的 BAM 文件是以 ${i}.bam 的形式存儲(chǔ)在 02Alignment 目錄中。
        #-o 選項(xiàng)指定了排序后的 BAM 文件的輸出路徑和文件名。排序后的 BAM 文件將保存為 ./02Alignment/${i}_sorted.bam。${i} 變量也在輸出文件名中使用，以保持結(jié)果文件的命名與輸入文件相關(guān)。
    samtools sort -@ ${threads} ./02Alignment/${i}.bam -o ./02Alignment/${i}_sorted.bam  
    
    #為了節(jié)約空間，可以將不用的bam文件刪除
    sudo rm -r ./02Alignment/${i}.bam
    
    #這是 stringtie 工具的名稱，它是一個(gè)用于轉(zhuǎn)錄組裝和量化的常用工具。
    #-p 選項(xiàng)指定了要使用的線程數(shù)，這里設(shè)置為 22，表示使用 22 個(gè)線程來(lái)執(zhí)行裝配和量化任務(wù)，以提高處理速度。
    #-G 選項(xiàng)指定了參考基因組的注釋文件的路徑，這是一個(gè) GTF 文件。stringtie 使用這個(gè)文件來(lái)輔助裝配和注釋轉(zhuǎn)錄本。
    #-e 選項(xiàng)表示要生成輸入 BAM 文件中每個(gè)樣本的單個(gè)轉(zhuǎn)錄本集合。
    #-B 選項(xiàng)表示生成一個(gè)包含每個(gè)樣本的 Ballgown 格式輸出。Ballgown 是一種用于分析和可視化 RNA-seq 數(shù)據(jù)的工具，可以用于分析基因表達(dá)和調(diào)查轉(zhuǎn)錄本的變化。
    #-o 選項(xiàng)指定了裝配和量化后的轉(zhuǎn)錄本 GTF 文件的輸出路徑和文件名。${i} 變量似乎用于迭代多個(gè)樣本文件，因此 ${i} 會(huì)在循環(huán)中被實(shí)際的文件名替換。這里生成的 GTF 文件將保存在 
    # -A 選項(xiàng)用于生成基因豐度的 TSV 文件，記錄每個(gè)樣本中每個(gè)基因的豐度信息。這個(gè)文件的路徑和文件名也包括 ${i} 變量，以保持結(jié)果文件的命名與輸入文件相關(guān)。生成的 TSV 文件將保存在 ./ballgown/${i}/gene_abundances.tsv。
    # 這是輸入 BAM 文件的路徑，即已經(jīng)比對(duì)和排序的測(cè)序數(shù)據(jù)的文件。${i} 變量也在輸入文件名中使用，以處理多個(gè)樣本文件。
    stringtie -p ${threads} -G ${GTF} -e -B -o ./ballgown/${i}/transcripts.gtf -A ./ballgown/${i}/gene_abundances.tsv ./02Alignment/${i}_sorted.bam 

done    
#計(jì)算reads數(shù)
featureCounts -T ${threads} -p -a ${GTF} -o ./FeatureCounts_gene.txt -t gene -g gene_id  ./02Alignment/*_sorted.bam 


# 提取stringtie產(chǎn)生的FPKM以及TPM值
Rscript -e "rm(list = ls()); setwd('./'); 
library(dplyr); library(stringr); library(data.table);

myresult = list();

myresult$FPKM = list.files('ballgown/') %>% lapply(function(x){ 
  file = read.table(paste0('ballgown/', x, '/gene_abundances.tsv'), sep = '\t', header = T, quote = '', check.names = F) %>% select(c('Gene ID', 'FPKM'))
  colnames(file)[2] = x; return(file)
}) %>% Reduce(x = ., function(x, y){ full_join(x, y, by = 'Gene ID') });

myresult$TPM = list.files('ballgown/') %>% lapply(function(x){ 
  file = read.table(paste0('ballgown/', x, '/gene_abundances.tsv'), sep = '\t', header = T, quote = '', check.names = F) %>% select(c('Gene ID', 'TPM'))
  colnames(file)[2] = x; return(file)
}) %>% Reduce(x = ., function(x, y){ full_join(x, y, by = 'Gene ID') });

openxlsx::write.xlsx(myresult, 'Counts.tpm_fpkm.xlsx', rowNames = F, overwrite = T);
"