劉小澤寫于18.8.21
爭取用三部分整理完
Part I是前情提要齐苛，學習一些背景知識炼团；
Part II是實戰(zhàn)，從文章拿到數(shù)據(jù)開始承二，只要是服務器能完成的工作都完成；
Part III 將是R處理部分纲爸，也是最核心的部分亥鸠，我會先學習Bioconductor，抓緊時間吧
. 閱讀文獻
文章梗概
運動神經(jīng)元存活蛋白(Survival Motor Neuron, SMN)是人體內(nèi)的泛表達蛋白识啦，它的會導致脊髓性肌萎速癥(Spinal muscular atrophy, SMA)负蚊，SMA是由SMN1基因突變引起的，但是這個基因是廣泛表達的基因颓哮，為什么運動神經(jīng)元(Motor Neurons, MNs)偏偏是最易受影響的細胞類型之一家妆？實驗通過對照組和SMA患者誘導多能干細胞（IPSC）從而產(chǎn)生純化的運動神經(jīng)元，通過固定冕茅、抗體標記伤极，然后進行了RNA測序研究。在患者運動神經(jīng)元中發(fā)現(xiàn)了SMA特異性的變化姨伤，其中包括內(nèi)質網(wǎng)(ER)應激通路的過度激活哨坪。功能研究表明，抑制內(nèi)質網(wǎng)的應激反應能提高運動神經(jīng)元的存活率乍楚。在患病小鼠中当编，使用ER應激抑制劑穿過血腦屏障可以保護脊髓運動神經(jīng)元。因此徒溪，選擇性激活內(nèi)質網(wǎng)應激反應忿偷，導致了SMA患病個體的運動神經(jīng)元死亡。
開始實戰(zhàn)
萬事之源臊泌，軟件為先
#配置conda
conda create -n rna-seq python=2 samtools fastp sra-tools hisat2 rseqc -y
conda install -c bcbio htseq -y
conda install numpy pysam -y
CONBIN=/home/biosoft/miniconda3/envs/rna-seq/bin
配置好工作路徑
RNA=/home/project/rna-seq
mkdir -p $RNA/{raw_data,clean_data,ref/{genome,gtf,index},align,stats,matrix}
RAW=$RNA/raw_data
CLEAN=$RNA/clean-data
GENOME=$RNA/ref/genome
GTF=$RNA/ref/gtf
INDEX=$RNA/ref/index
ALIGN=$RNA/aign
MATRIX=$RNA/matrix
STATS=$RNA/stats
下載測試數(shù)據(jù)
GEO數(shù)據(jù)庫中的編號是：GSE69175
關于數(shù)據(jù)：
NGS測序數(shù)據(jù)一般會上傳到SRA數(shù)據(jù)庫里面鲤桥，但是怎么從GEO數(shù)據(jù)庫中找到SRA原始數(shù)據(jù)的下載地址？【GEO數(shù)據(jù)庫地址：https://www.ncbi.nlm.nih.gov/geo/】
具體的層級關系是：SRP(項目)—>SRS(樣本)—>SRX(數(shù)據(jù)產(chǎn)生)—>SRR(數(shù)據(jù)本身)
SRR數(shù)據(jù)庫地址：https://ftp-private.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/
##################################################################
################# Part 1: Data download #########################
##################################################################
# SRR2038215-SRR2038216: 對照
# SRR2038217-SRR2038218: SMN
cd $RAW
for i in `seq 15 18`;do
ascp -v -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
-k 1 -T -l200m anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR358/SRR35899${i}/SRR35899${i}.sra ./ #使用ascp下載缺虐，嗖嗖嗖～
done
fastq-dump --gzip --split-3 *.sra #sra轉為fastq
#以下3行原版來自生信技能樹Jimmy
find $RAW -name "*.gz" | grep 1.fastq.gz >1
find $RAW -name "*.gz" | grep 2.fastq.gz >2
paste 1 2 > raw_conf #將read1和read2各自的合集再整合到一起芜壁，形成一個數(shù)據(jù)配置文件
cp raw_conf $CLEAN
質控過濾
這里使用的是fastp，同時融合了質控、過濾的模式慧妄；
同時也可以使用fastqc + trimmomatic/ trim_galore進行
##################################################################
################## Part 2: QC & trim ############################
##################################################################
source activate rna-seq
cd $CLEAN
#下面四行原版來自生信技能樹Jimmy
cat raw_conf | while read id;do
line=(${id})
fq1=${line[0]}
fq2=${line[1]}
fastp -i $fq1 -I $fq2 -o out.$(basename $fq1) -O out.$(basename $fq2)-w 16
done
source deactivate
下載參考基因組和注釋文件
人類參考基因組選擇UCSC數(shù)據(jù)庫顷牌；注釋文件選擇GENCODE,https://www.gencodegenes.org/
下載好基因組，需要構建基因組索引
如果自己的項目做非模式物種塞淹，可以用二代三代測序數(shù)據(jù)組裝成參考轉錄組窟蓝，例如trinity組裝。如果做的物種有參考基因組饱普，就方便一些运挫，可以直接從相關數(shù)據(jù)庫中下載參考基因組
##################################################################
############### Part 3: prepare ref & index ######################
##################################################################
##Download hg19
cd $GENOME
for i in $(seq 1 22) X Y M;do
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz;
done
gunzip *.gz
for i in $(seq 1 22) X Y M;do
cat chr${i}.fa >> hg19.fa
done
rm chr*
##或者也可以直接下載成品
#wget http://10.10.0.8/hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz
##Download GRCh38 https://www.gencodegenes.org/releases/current.html
cd $GTF
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/GRCh37_mapping/gencode.v28lift37.annotation.gtf.gz
gunzip *.gz
## hisat2 index
cd $INDEX
hisat2-build -p 20 $GENOME/hg19.fa hg19 #-p為線程數(shù)
比對
##################################################################
######## Part 4: align & sort & index human samples:56-58 ########
##################################################################
source activate rna-seq
cd $ALIGN
for i in `seq 15 18`;do
hisat2 -t -p 20 -x $INDEX/hg19 \
-1 $CLEAN/out.SRR20382${i}_1.fastq.gz \
-2 $CLEAN/out.SRR20382${i}_2.fastq.gz \
-S SRR20382${i}.sam
samtools view -Sb SRR20382${i}.sam > SRR20382${i}.bam
samtools sort -@ 20 -o SRR20382${i}.sorted.bam SRR20382${i}.bam
samtools index SRR20382${i}.sorted.bam
rm *.sam
done
source deactivate
#關于排序：默認按染色體位置排序；-n根據(jù)read名排序套耕；-t根據(jù)tag排序
使用了samtools的三件套：轉換(view)谁帕、排序(sort)、建索引(index)
轉換: -S指輸入文件格式（不加-S默認輸入是bam）冯袍，-b指定輸出文件（默認輸出sam）【如果要bam轉sam匈挖，-h設置輸出sam時帶上頭注釋信息】
排序： 對bam排序，-@ 線程數(shù)康愤， -o輸出文件
索引： 結果會產(chǎn)生.bai文件【必須排序后才能建索引～就像體育課必須從高到矮排好以后再報數(shù)】
基本信息統(tǒng)計
##################################################################
################ Part 5: basic statistics #######################
##################################################################
source activate rna-seq
cd $STATS
for i in `seq 15 18`;do
samtools flagstat $ALIGN/SRR20382${i}.sorted.bam > SRR20382${i}.sorted.flag
done
#如果想根據(jù)flag提取特定區(qū)域儡循，比如想查看1號染色體100-10000的區(qū)域的信息
#samtools view -b -f 0x10 $ALIGN/SRR20382${i}.sorted.bam chr1:100-10000 > ${i}.flag.bam
#samtools flagstat ${i}.flag.bam
#使用RSeqQC統(tǒng)計
#先用bam_stat.py對bam文件統(tǒng)計，看下比對率
#bam_stat.py -i $ALIGN/SRR20382${i}.sorted.bam > SRR20382${i}.bam.stat
source deactivate
reads計數(shù)
基于基因組水平征冷，可以用Htseq-count和featureCounts
-
Htseq-count：它是用python寫的一個腳本择膝，
conda install -c bcbio htseq -y安裝好以后可以直接拿來用需要比對文件sam/bam格式；需要基因組注釋文件gff/gtf格式检激。這兩個文件的染色體名稱必須相同肴捉，因為需要將比對位置和特征信息相匹配就要借助染色體名。
它的工作原理是：先通過bam文件找到reads比對上的外顯子呵扛，然后去gtf文件中將外顯子的基因ID進行統(tǒng)計每庆，得到的結果就是未經(jīng)標準化的表達量source activate rna-seq cd $MATRIX for i in `seq 15 18`;do $CONBIN/htseq-count -s no -r name -f bam $ALIGN/SRR20382${i}.sorted.bam \ $GTF/gencode.v28lift37.annotation.gtf \ >SRR20382${i}.count 2>SRR20382${i}.log done source deactivatehtseq-count的
-s參數(shù)是strand意思，默認使用鏈特異性建庫方式今穿，也就是計算過程中要求read1只能和雙端測序的其中一條比較，read2只能和另一條比較伦籍；我們一般設置no即可蓝晒，表示每條read都和正負鏈比較一次；
-r參數(shù)兩個選擇name和pos帖鸦，數(shù)據(jù)根據(jù)位置或者read名稱排序芝薇，默認name；
-f參數(shù)輸入文件格式bam/sam作儿；-m參數(shù)一般也就是默認union模式結果每個樣本輸出一個count文件洛二，其中包含了基因名和reads計數(shù)；另外，如果你看到文件倒數(shù)幾行晾嘶，會發(fā)現(xiàn)還有幾行帶文字的
htseq-count計數(shù)no_feature：比對區(qū)域與任何基因都沒有重疊
ambiguous：比對區(qū)域與多個基因都發(fā)生重疊
too_low_aQual：比對質量低于設定閾值（默認是10）
not_aligned：無法比對上
alignment_not_unique：比對位置不唯一 -
featureCounts：隸屬于subread套件【相比于htseq更快】
source activate rna-seq cd $MATRIX begin=$(date +%s) for i in `seq 15 18`;do featureCounts -T 20 -p -t exon -g gene_id -a $GTF/gencode.v28lift37.annotation.gtf -o SRR20382${i}.feature.count $ALIGN/SRR20382${i}.bam done tim=$(echo "$(date +%s)-$begin" | bc) printf "time of featureCounts for 4 samples: %.4f seconds" $tim source deactivate #運行結果顯示： #time of featureCounts for 4 samples: 366.5310 seconds-T參數(shù)設置線程妓雾，
-p參數(shù)表示針對雙端測序數(shù)據(jù)，
-t 參數(shù)默認將外顯子認定為基因組的一個feature(搜索特征)垒迂，
-g參數(shù)默認按照基因名來計數(shù)械姻，
-a參數(shù)輸入注釋文件，
-o參數(shù)輸出文件結果有兩個机断，一個count文件楷拳，一個summary文件
featureCounts文件
對兩個軟件的結果進行合并
##合并htdeq生成的count文件到matrix.count
cd $MATRIX/htseq
perl -lne 'if ($ARGV=~/(.*).count/){print "$1\t$_"}' *.count | grep -v "_" >matrix.count
##合并featureCounts生成的count文件到matrix_2.count
cd $MATRIX/featureCounts
for i in `seq 15 18`;do
cut -f 1,7 SRR20382${i}.feature.count | grep -v "^#" > SRR20382${i}.matrix
sed -i '1d' SRR20382${i}.matrix
done
perl -lne 'if ($ARGV=~/(.*).matrix/){print "$1\t$_"}' *.matrix >matrix_2.count
統(tǒng)計一下兩個軟件的計數(shù)之和
#統(tǒng)計featureCounts
perl -alne '$sum += $F[2]; END {print $sum}' matrix.count
#結果是5882943
#統(tǒng)計htseq-count,結果是786338
歡迎關注我們的公眾號～_～
我們是兩個農(nóng)轉生信的小碩，打造生信星球吏奸，想讓它成為一個不拽術語欢揖、通俗易懂的生信知識平臺。需要幫助或提出意見請后臺留言或發(fā)送郵件到Bioplanet520@outlook.com
