暴躁版
制約我們的是熱情嗎旧巾？是知識嗎耸序？
不，是網(wǎng)速啊鲁猩，阿西吧坎怪！
1 獲取SRA號——自己讀文獻(xiàn)，去pubmend搜
2 下載數(shù)據(jù)
cat SRR_Acc_List.txt | while read id; do (prefetch ${id} );done
別試了廓握，龜速搅窿，下不下來，請愉快的打上一個&然后洗洗睡吧
jimmy貼心的給了一個kill命令的代碼
ps -ef | grep prefetch | awk '{print $2}' | while read id; do kill ${id}; done
呵呵噠
3 質(zhì)控
sra轉(zhuǎn)fastq
for i in $wkd/*sra
do
echo $i
fastq-dump --split-3 --skip-technical --clip --gzip $i ## 批量轉(zhuǎn)換
done
核心：fastq-dump
--split-3 確保一個雙端測序的樣本被拆分成兩個fastq文件隙券。
--skip-technical Dump only biological reads
-W|--clip Remove adapter sequences from reads
--gzip Compress output using gzip: deprecated, not
recommended （哈哈男应，為啥不推薦呢）
fastq轉(zhuǎn)fastqc
ls *gz | xargs fastqc -t 4
multiqc ./
-t 線程數(shù) （double哈哈）
還是好慢，rna-seq的正確打開方法是娱仔，找本書沐飘，一邊seq一邊看
multiqc一下，就得到了multiqc的報(bào)告，下下來看一下耐朴，打不開借卧，loading report ，F(xiàn)U...K!馬景濤咆哮版筛峭！
4.過濾
過濾質(zhì)量差的reads和去除接頭
ls ~/sra/*_1.fastq.gz >1
ls ~/sra/*_2.fastq.gz >2
paste 1 2 >config
cat 一下config,長這樣
/trainee2/hz10/sra/SRR1039510_1.fastq.gz /trainee2/hz10/sra/SRR1039510_2.fastq.gz
/trainee2/hz10/sra/SRR1039511_1.fastq.gz /trainee2/hz10/sra/SRR1039511_2.fastq.gz
/trainee2/hz10/sra/SRR1039512_1.fastq.gz /trainee2/hz10/sra/SRR1039512_2.fastq.gz
trim_galore ,用于去除低質(zhì)量和接頭數(shù)據(jù)
改一下大神的腳本铐刘，換成自己的路徑，長這樣
conda activate rna2
bin_trim_galore=trim_galore
dir='/trainee2/hz10/sra/clean'
cat $1 |while read id
do
arr=(${id})
fq1=${arr[0]}
fq2=${arr[1]}
$bin_trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o $dir $fq1 $fq2
done
conda deactivate rna2
然后bash一下蜒滩，好慢滨达，我忘了打&，然后掉線了俯艰，我恨捡遍！我需要把之前trim好的刪掉再重新bash嗎？還是它會很智能的自己跳過去呢竹握？它好像跳過去了画株，nice~
太慢了，回家吃飯了啦辐。谓传。。芹关。续挟。
好了，如下侥衬。
├── [ 123] 1
├── [ 123] 2
├── [ 246] config
├── [ 284] qc.sh
├── [2.9K] SRR1039510_1.fastq.gz_trimming_report.txt
├── [1.2G] SRR1039510_1_val_1.fq.gz
├── [3.1K] SRR1039510_2.fastq.gz_trimming_report.txt
├── [1.2G] SRR1039510_2_val_2.fq.gz
├── [2.9K] SRR1039511_1.fastq.gz_trimming_report.txt
├── [1.2G] SRR1039511_1_val_1.fq.gz
├── [3.1K] SRR1039511_2.fastq.gz_trimming_report.txt
├── [1.2G] SRR1039511_2_val_2.fq.gz
├── [2.9K] SRR1039512_1.fastq.gz_trimming_report.txt
├── [1.5G] SRR1039512_1_val_1.fq.gz
├── [3.1K] SRR1039512_2.fastq.gz_trimming_report.txt
└── [1.5G] SRR1039512_2_val_2.fq.gz
--phred33 Instructs Cutadapt to use ASCII+33 quality scores as Phred scores (Sanger/Illumina 1.9+ encoding) for quality trimming.
--stringency Overlap with adapter sequence required to trim a sequence. 剪切掉的overlap的序列
***--length <INT> **** Discard reads that became shorter than length INT because of either quality or adapter trimming. A value of '0' effectively disables this behaviour. Default: 20 bp.
--paired 雙端測序 并且需要雙端測序的兩條序列都比length的設(shè)定值要長诗祸，可以再不影響fastq格式的情況下，過濾掉過短的序列轴总。
4. 比對
使用hisat2
構(gòu)建hisat2基因組索引,來自生信技能樹論壇直颅，傳說很慢
# 構(gòu)建目錄
mkdir /mnt/d/reference/index
mkdir /mnt/d/reference/index/hisat
#下載
wget -P /mnt/d/reference/index/hisat [url=ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz]ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz[/url]
#解壓，-C指定解壓目錄
tar -zxvf /mnt/d/reference/index/hisat/hg19.tar.gz -C /mnt/d/reference/index/hisat
# 移除下載的壓縮包
rm /mnt/d/reference/index/hisat/*.tar.gz
#查看解壓文件
ll /mnt/d/reference/index/hisat/hg19
Tips_1:不用單獨(dú)給下載的hg19的index文件再分別構(gòu)建目錄怀樟，解壓的時候會自動創(chuàng)建hg19目錄功偿。
Tips_2: 解壓的文件中，包含genome.*.ht2的8個文件和一個shell腳本往堡。
hisat2用法
hisat2 [options]* -x <hisat2-idx>{-1 <m1> -2 <m2> | -U <r> } [-S <hit>]
-x 指定基因組索引
-1 指定第一個fastq文件:雙端測序結(jié)果的第一個文件械荷。若有多組數(shù)據(jù)，使用逗號將文件分隔虑灰。Reads的長度可以不一致养葵。
-2 指定第二個fastq文件:雙端測序結(jié)果的第二個文件。若有多組數(shù)據(jù)瘩缆，使用逗號將文件分隔，并且文件順序要和-1參數(shù)對應(yīng)佃蚜。Reads的長度可以不一致庸娱。
-S 指定輸出的SAM文件着绊。
改了老師的批量運(yùn)行代碼
cd $wkd/clean
ls *gz|cut -d"_" -f 1 |sort -u |while read id;do
ls -lh ${id}_1_val_1.fq.gz ${id}_2_val_2.fq.gz
hisat2 -p 10 -x ~/sra/index/hisat/hg38/genome -1 ${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz -S ${id}.hisat.sam
done
運(yùn)行了一下，運(yùn)行完第一個不錯熟尉，比對98%归露，第二個卡住了，報(bào)錯斤儿，google了一下原因剧包，sam文件太大，內(nèi)存不夠往果，服務(wù)器爆掉了疆液。。陕贮。堕油。我太難了。肮之。掉缺。。
然后就一個一個做戈擒，運(yùn)行完了轉(zhuǎn)sam眶明，或者直接轉(zhuǎn)sam。
hisat2 -p 1 -x ~/sra/index/hisat2/hg38/genome -1 SRR1039511_1_val_1.fq.gz -2 SRR1039511_2_val_2.fq.gz | samtools sort -@ 1 -o ~/sra/clean/SRR1039511.sort.bam -
bam 轉(zhuǎn) sam
samtools sort -O bam -@ 5 -o SRR1039510.hisat.bam SRR1039510.hisat.sam
-O 輸出文件格式
-@ 線程數(shù)
最后是要操作的文件名筐高，上面用管道符的是最后的-
我想買個服務(wù)器搜囱，我太難了。凯傲。犬辰。。
終于結(jié)束了冰单，為bam文件建立索引
ls *.bam |xargs -i samtools index {}
ls *.bam |while read id ;do ( samtools flagstat -@ 1 $id > $(basename ${id} ".bam").flagstat );done
samtools view SAM轉(zhuǎn)換為二進(jìn)制對應(yīng)的BAM格式
samtools sort 比對排序
samtools index 構(gòu)建索引幌缝。對排序文件進(jìn)行索引之后，有利于快速提取基因組重疊區(qū)域的比對結(jié)果
samtools flagstats 給出BAM文件的比對結(jié)果
samtools常用命令詳解
SRR1039512.sort.bam.bam.flagstat文件可以直接cat查看比對結(jié)果
5.獲得表達(dá)矩陣
gtf="/teach/database/gtf/gencode.v29.annotation.gtf.gz"
featureCounts -T 5 -p -t exon -g gene_id -a $gtf -o all.id.txt ~/sra/clean/*.bam 1>counts.id.log 2>&1 &
featurecounts 據(jù)說最大的優(yōu)點(diǎn)是快
-a 輸入GTF/GFF基因組注釋文件
-p 這個參數(shù)是針對paired-end數(shù)據(jù)诫欠。Check validity of paired-end distance when counting read pairs. Use -d and -D to set thresholds.
-F 指定-a注釋文件的格式涵卵，默認(rèn)是GTF
-g 從注釋文件中提取Meta-features信息用于read count，默認(rèn)是gene_id
-t 跟-g一樣的意思荒叼，其是默認(rèn)將exon作為一個feature
-o 輸出文件
-T 多線程數(shù)
然后我們就得到表達(dá)矩陣 all.id.txt 啦轿偎，就可以用R進(jìn)行下游分析啦啦啦啦（啦啦啦你妹），limma,DESeq2被廓，想用什么用什么坏晦，哪里不會點(diǎn)哪里~
ps，我覺得jimmy大佬有一個地方寫錯了，align下面沒有bam文件昆婿，應(yīng)該是$wkd/clean/*bam~ 覺得自己超棒球碉，有沒有~
一天的時間，跑完了一個流程仓蛆，但是睁冬，只跑了3個樣本，服務(wù)器就over了看疙，sad
明天再跟著洲更大神的視頻跑一個流程豆拨，我是不是就算入門了呢？
（不交錢是我最后的倔強(qiáng)~）
附：jimmy原教程
