1 ###文件下載
單個(gè)文件
wget -c -t 0 -O SRR13107018.sra https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR13107018/SRR13107018
#-c -t 配合使用可以防止下載數(shù)據(jù)的過(guò)程中鏈接中斷的問(wèn)題捕仔，-O則可以指定下載路徑和文件名。
#多個(gè)文件一起下
/home/software/sratoolkit.2.9.6-1-centos_linux64/bin/prefetch -O ./ --option-file SRR_Acc_List.txt
2 ## 循環(huán)把下載的所有sra文件都變?yōu)閒astq(雙端測(cè)序數(shù)據(jù))
ls *sra | while read id ;do (nohup /home/software/sratoolkit.2.9.6-1-centos_linux64/bin/fastq-dump --gzip --split-3 $id &) ;done
3 ### fastp質(zhì)量控制
ls *fastq.gz | cut -d"_" -f 1 |sort -u | while read id ;do (nohup fastp -i ${id}_1.fastq.gz -I ${id}_2.fastq.gz -g -q 5 -n 15 -l 150 -u 50 -o ${id}_1.filter.fastq.gz -O ${id}_2.filter.fastq.gz -h ${id}.fastp.html &) ;done
# 去除帶接頭(adapter)的雙端讀長(zhǎng)(默認(rèn)開啟风宁，可用-A關(guān)閉)；
# 去掉單端讀長(zhǎng)中含 N(N 表示無(wú)法確定堿基信息)的堿基比例大于10%部分（默認(rèn)開啟）
# -q 設(shè)置堿基質(zhì)量閾值曙旭，默認(rèn)是15茧球，phred質(zhì)量評(píng)分≥Q15
# -n 一條read中N堿基的數(shù)量超過(guò)了多少個(gè)，那么這條read就被刪除枪狂，默認(rèn)是5（即這條read里有5個(gè)N）
# -g 進(jìn)行PolyG尾的去除
# -u 允許百分之多少的堿基不合格（0-100）危喉，默認(rèn)是40（就是說(shuō)40%），超過(guò)這個(gè)比例州疾，整條read就被刪除
# -l read小于這個(gè)參數(shù)設(shè)定長(zhǎng)度會(huì)被丟棄或刪除辜限，默認(rèn)是15,由你的測(cè)序bp決定
# -j, --json 輸出的json報(bào)告文件名，以“.json”結(jié)尾
# -h, --html 輸出的html報(bào)告文件名严蓖，以“.html”結(jié)尾
4 ###參考基因組下載及建索引（如已有薄嫡，忽略此步）
## 下載基因組序列
mkdir -p genome/ARS-UCD1.2 && cd genome/ARS-UCD1.2
nohup wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/263/795/GCF_002263795.1_ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.fna.gz &
gunzip GCF_002263795.1_ARS-UCD1.2_genomic.fna.gz
## 下載gtf文件并解壓
wget https://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_mammalian/Bos_taurus/latest_assembly_versions/GCF_002263795.1_ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.gtf.gz
gunzip GCF_002263795.1_ARS-UCD1.2_genomic.gtf.gz
## 建索引（hisat2軟件）
mkdir index & cd index
nohup hisat2-build -p 4 GCF_002263795.1_ARS-UCD1.2_genomic.fna GCF_002263795.1_ARS-UCD1.2_genomic > hisat2.log &
5 ###reads mapping到參考基因組-----hisat2
mkdir hismap.sam
ls *filter.fastq.gz|cut -d"_" -f 1 |sort -u | while read id ;do (nohup hisat2 -p 8 -x /home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic -1 ${id}_1.filter.fastq.gz -2 ${id}_2.filter.fastq.gz -S hismap.sam/${id}.hismap.sam &) ;done
#綿羊參考基因組
/home/sll/genome-sheep/Oar_rambouillet_v1.0-ncbi/GCF_002742125.1_Oar_rambouillet_v1.0_genomic
#牛參考基因組
/home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic
# -p: 線程數(shù)
# -x: 基因組索引前綴氧急。所下的基因組索引為多個(gè)文件，索引前綴到genome為止毫深。
# -1/-2: fastq輸入文件吩坝。當(dāng)輸入為單端測(cè)序時(shí)使用-U 指定輸入。單端或雙端的輸入都可使用逗號(hào)分隔輸入多個(gè)樣本哑蔫，或使用多次-1 -2 / -U 指定多個(gè)輸入钉寝。e.g.: -U sample1.fq.gz,sample2.fq.gz -U sample3.fq.gz
# -S: 輸出sam文件路徑。
6 ###轉(zhuǎn)為bam文件并排序
cd hisout.sam
ls *sam | cut -d"." -f 1 |sort -u | while read id ;do (nohup samtools view -S ${id}.hismap.sam -b > ${id}.hismap.bam &) ;done
ls *hismap.bam | cut -d"." -f 1 |sort -u | while read id ;do (nohup samtools sort -@ 8 ${id}.hismap.bam -o ${id}_sort.bam &) ;done
# sort: 進(jìn)行排序
# -@: 線程數(shù)
# -o: 輸出bam文件
# 最后一項(xiàng)為輸入文件
7 ###為sort.bam文件建索引
ls *sort.bam | cut -d"_" -f 1 | sort -u | while read id ;do (nohup samtools index ${id}_sort.bam ${id}_sort.bam.index &) ;done
## index用于建立索引闸迷，使之快速訪問(wèn)bam文件
8 ###faturecount計(jì)數(shù)定量
nohup featureCounts -p -t exon -g gene_id -M -T 8 -a /home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.gtf -o all.featurecounts.txt *_sort.bam &
#綿羊gtf文件
/home/sll/genome-sheep/Oar_rambouillet_v1.0-ncbi/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.gtf
#牛gtf文件
/home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.gtf
# featureCounts 進(jìn)行定量嵌纲，統(tǒng)計(jì)比對(duì)在這個(gè)基因的坐標(biāo)上的read數(shù)
# -t 設(shè)置feature-type，-t指定的必須是gtf中有的feature稿黍，同時(shí)read只有落到這些feature上才會(huì)被統(tǒng)計(jì)到疹瘦，默認(rèn)是“exon”
# -p 雙端數(shù)據(jù)時(shí)需要
# -a 轉(zhuǎn)錄組注釋文件
# -o 輸出文件
###報(bào)錯(cuò)處理
cat GCF_002263795.1_ARS-UCD1.2_genomic.gtf | grep gene_id | wc -w #檢查原始gtf注釋文件gene_id個(gè)數(shù)
cat GCF_002263795.1_ARS-UCD1.2_genomic.gtf | grep gene_id\ \"\"\; | wc -w #檢查空gene_id值的行數(shù)
sed -i -e '/gene_id\ \"\"\;/d' GCF_002263795.1_ARS-UCD1.2_genomic.gtf #刪除包括空gene_id值的行
cat GCF_002263795.1_ARS-UCD1.2_genomic.gtf | grep gene_id\ \"\"\; | wc -w #修改后的檢查一下 檢查包括空gene_id值的行 0為已刪除掉
9 ###對(duì)featureCounts的結(jié)果進(jìn)行整合，html文件可視化
multiqc all.featurecounts.txt.summary -o all.counts.summary
10 ###提取定量信息
awk -F '\t' '{print $1,$7,$8,$9,$10,$11,$12}' OFS='\t' all.featurecounts.txt > all_fcount.matrix.txt
# 1列為基因id巡球，7列及以后為樣本定量信息
# \t 表示以制表符分割開來(lái)
11 ###將矩陣導(dǎo)入R中言沐，采用DESeq2進(jìn)行差異分析
